diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index 11e464ebe..fbd84f695 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,1548 +1,1556 @@
/**
* @file aka_types.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 17 2011
* @date last modification: Tue Feb 20 2018
*
* @brief description of the "simple" types
*
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu 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.
*
* Akantu 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 Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_math.hh"
/* -------------------------------------------------------------------------- */
#include <initializer_list>
#include <iomanip>
#include <type_traits>
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_AKA_TYPES_HH_
#define AKANTU_AKA_TYPES_HH_
namespace akantu {
enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) };
/**
* DimHelper is a class to generalize the setup of a dim array from 3
* values. This gives a common interface in the TensorStorage class
* independently of its derived inheritance (Vector, Matrix, Tensor3)
* @tparam dim
*/
template <UInt dim> struct DimHelper {
static inline void setDims(UInt m, UInt n, UInt p,
- std::array<UInt, dim> dims);
+ std::array<UInt, dim> & dims);
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<1> {
static inline void setDims(UInt m, UInt /*n*/, UInt /*p*/,
- std::array<UInt, 1> dims) {
+ std::array<UInt, 1> & dims) {
dims[0] = m;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<2> {
static inline void setDims(UInt m, UInt n, UInt /*p*/,
- std::array<UInt, 2> dims) {
+ std::array<UInt, 2> & dims) {
dims[0] = m;
dims[1] = n;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<3> {
- static inline void setDims(UInt m, UInt n, UInt p, std::array<UInt, 3> dims) {
+ static inline void setDims(UInt m, UInt n, UInt p,
+ std::array<UInt, 3> & dims) {
dims[0] = m;
dims[1] = n;
dims[2] = p;
}
};
/* -------------------------------------------------------------------------- */
template <typename T, UInt ndim, class RetType> class TensorStorage;
/* -------------------------------------------------------------------------- */
/* Proxy classes */
/* -------------------------------------------------------------------------- */
namespace tensors {
template <class A, class B> struct is_copyable {
enum : bool { value = false };
};
template <class A> struct is_copyable<A, A> {
enum : bool { value = true };
};
template <class A> struct is_copyable<A, typename A::RetType> {
enum : bool { value = true };
};
template <class A> struct is_copyable<A, typename A::RetType::proxy> {
enum : bool { value = true };
};
} // namespace tensors
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
namespace types {
namespace details {
template <typename reference_> class vector_iterator {
public:
using difference_type = std::ptrdiff_t;
using value_type = std::decay_t<reference_>;
using pointer = value_type *;
using reference = reference_;
using iterator_category = std::input_iterator_tag;
vector_iterator(pointer ptr) : ptr(ptr) {}
// input iterator ++it
vector_iterator & operator++() {
++ptr;
return *this;
}
// input iterator it++
vector_iterator operator++(int) {
auto cpy = *this;
++ptr;
return cpy;
}
vector_iterator & operator+=(int n) {
ptr += n;
return *this;
}
vector_iterator operator+(int n) {
vector_iterator cpy(*this);
cpy += n;
return cpy;
}
// input iterator it != other_it
bool operator!=(const vector_iterator & other) const {
return ptr != other.ptr;
}
bool operator==(const vector_iterator & other) const {
return ptr == other.ptr;
}
difference_type operator-(const vector_iterator & other) const {
return this->ptr - other.ptr;
}
// input iterator dereference *it
reference operator*() { return *ptr; }
pointer operator->() { return ptr; }
private:
pointer ptr;
};
} // namespace details
} // namespace types
/**
* @class TensorProxy aka_types.hh
* The TensorProxy class is a proxy class to the
* TensorStorage it handles the wrapped case. That is to say if an accessor
* should give access to a Tensor wrapped on some data, like the
* Array<T>::iterator they can return a TensorProxy that will be automatically
* transformed as a TensorStorage wrapped on the same data
* @tparam T stored type
* @tparam ndim order of the tensor
* @tparam _RetType real derived type
*/
template <typename T, UInt ndim, class RetType_>
class TensorProxy : public TensorProxyTrait {
protected:
using RetTypeProxy = typename RetType_::proxy;
constexpr TensorProxy(T * data, UInt m, UInt n, UInt p) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->values = data;
}
template <class Other, typename = std::enable_if_t<
tensors::is_copyable<TensorProxy, Other>::value>>
explicit TensorProxy(const Other & other) {
this->values = other.storage();
for (UInt i = 0; i < ndim; ++i) {
this->n[i] = other.size(i);
}
}
public:
using RetType = RetType_;
UInt size(UInt i) const {
AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
<< " dimensions, not "
<< (i + 1));
return n[i];
}
inline UInt size() const {
UInt _size = 1;
for (UInt d = 0; d < ndim; ++d) {
_size *= this->n[d];
}
return _size;
}
T * storage() const { return values; }
template <class Other, typename = std::enable_if_t<
tensors::is_copyable<TensorProxy, Other>::value>>
inline TensorProxy & operator=(const Other & other) {
AKANTU_DEBUG_ASSERT(
other.size() == this->size(),
"You are trying to copy two tensors with different sizes");
std::copy_n(other.storage(), this->size(), this->values);
return *this;
}
// template <class Other, typename = std::enable_if_t<
// tensors::is_copyable<TensorProxy, Other>::value>>
// inline TensorProxy & operator=(const Other && other) {
// AKANTU_DEBUG_ASSERT(
// other.size() == this->size(),
// "You are trying to copy two tensors with different sizes");
// memcpy(this->values, other.storage(), this->size() * sizeof(T));
// return *this;
// }
template <typename O> inline RetTypeProxy & operator*=(const O & o) {
RetType(*this) *= o;
return static_cast<RetTypeProxy &>(*this);
}
template <typename O> inline RetTypeProxy & operator/=(const O & o) {
RetType(*this) /= o;
return static_cast<RetTypeProxy &>(*this);
}
protected:
T * values;
std::array<UInt, ndim> n;
};
/* -------------------------------------------------------------------------- */
template <typename T> class VectorProxy : public TensorProxy<T, 1, Vector<T>> {
using parent = TensorProxy<T, 1, Vector<T>>;
using type = Vector<T>;
public:
constexpr VectorProxy(T * data, UInt n) : parent(data, n, 0, 0) {}
template <class Other> explicit VectorProxy(Other & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline VectorProxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
// inline VectorProxy<T> & operator=(const VectorProxy && other) {
// parent::operator=(other);
// return *this;
// }
using iterator = types::details::vector_iterator<T &>;
using const_iterator = types::details::vector_iterator<const T &>;
iterator begin() { return iterator(this->storage()); }
iterator end() { return iterator(this->storage() + this->size()); }
const_iterator begin() const { return const_iterator(this->storage()); }
const_iterator end() const {
return const_iterator(this->storage() + this->size());
}
/* ------------------------------------------------------------------------ */
T & operator()(UInt index) { return this->values[index]; };
const T & operator()(UInt index) const { return this->values[index]; };
};
template <typename T> class MatrixProxy : public TensorProxy<T, 2, Matrix<T>> {
using parent = TensorProxy<T, 2, Matrix<T>>;
using type = Matrix<T>;
public:
MatrixProxy(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
template <class Other> explicit MatrixProxy(Other & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline MatrixProxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
};
template <typename T>
class Tensor3Proxy : public TensorProxy<T, 3, Tensor3<T>> {
using parent = TensorProxy<T, 3, Tensor3<T>>;
using type = Tensor3<T>;
public:
Tensor3Proxy(const T * data, UInt m, UInt n, UInt k)
: parent(data, m, n, k) {}
Tensor3Proxy(const Tensor3Proxy & src) : parent(src) {}
Tensor3Proxy(const Tensor3<T> & src) : parent(src) {}
/* ---------------------------------------------------------------------- */
template <class Other>
inline Tensor3Proxy<T> & operator=(const Other & other) {
parent::operator=(other);
return *this;
}
};
/* -------------------------------------------------------------------------- */
/* Tensor base class */
/* -------------------------------------------------------------------------- */
template <typename T, UInt ndim, class RetType>
class TensorStorage : public TensorTrait {
public:
using value_type = T;
friend class Array<T>;
protected:
template <class TensorType> void copySize(const TensorType & src) {
for (UInt d = 0; d < ndim; ++d) {
this->n[d] = src.size(d); // NOLINT
}
this->_size = src.size();
}
TensorStorage() = default;
TensorStorage(const TensorProxy<T, ndim, RetType> & proxy) {
this->copySize(proxy);
this->values = proxy.storage();
this->wrapped = true;
}
public:
TensorStorage(const TensorStorage & src) = delete;
TensorStorage(const TensorStorage & src, bool deep_copy) : values(nullptr) {
if (deep_copy) {
this->deepCopy(src);
} else {
this->shallowCopy(src);
}
}
protected:
TensorStorage(UInt m, UInt n, UInt p, const T & def) {
static_assert(std::is_trivially_constructible<T>{},
"Cannot create a tensor on non trivial types");
DimHelper<ndim>::setDims(m, n, p, this->n);
this->computeSize();
this->values = new T[this->_size]; // NOLINT
this->set(def);
this->wrapped = false;
}
TensorStorage(T * data, UInt m, UInt n, UInt p) {
DimHelper<ndim>::setDims(m, n, p, this->n);
this->computeSize();
this->values = data;
this->wrapped = true;
}
public:
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void shallowCopy(const TensorType & src) {
this->copySize(src);
if (!this->wrapped) {
delete[] this->values;
}
this->values = src.storage();
this->wrapped = true;
}
/* ------------------------------------------------------------------------ */
template <class TensorType> inline void deepCopy(const TensorType & src) {
this->copySize(src);
if (!this->wrapped) {
delete[] this->values;
}
static_assert(std::is_trivially_constructible<T>{},
"Cannot create a tensor on non trivial types");
this->values = new T[this->_size]; // NOLINT
static_assert(std::is_trivially_copyable<T>{},
"Cannot copy a tensor on non trivial types");
std::copy_n(src.storage(), this->_size, this->values);
this->wrapped = false;
}
virtual ~TensorStorage() {
if (!this->wrapped) {
delete[] this->values;
}
}
/* ------------------------------------------------------------------------ */
inline TensorStorage & operator=(const TensorStorage & other) {
- if(this == &other) {
+ if (this == &other) {
return *this;
}
this->operator=(aka::as_type<RetType>(other));
return *this;
}
- inline TensorStorage & operator=(TensorStorage && other) noexcept = default;
+ // inline TensorStorage & operator=(TensorStorage && other) noexcept {
+ // std::swap(n, other.n);
+ // std::swap(_size, other._size);
+ // std::swap(values, other.values);
+ // std::swap(wrapped, other.wrapped);
+ // return *this;
+ // }
/* ------------------------------------------------------------------------ */
inline TensorStorage & operator=(const RetType & src) {
if (this != &src) {
if (this->wrapped) {
static_assert(std::is_trivially_copyable<T>{},
"Cannot copy a tensor on non trivial types");
// this test is not sufficient for Tensor of order higher than 1
AKANTU_DEBUG_ASSERT(this->_size == src.size(),
"Tensors of different size ("
<< this->_size << " != " << src.size() << ")");
std::copy_n(src.storage(), this->_size, this->values);
} else {
deepCopy(src);
}
}
return *this;
}
/* ------------------------------------------------------------------------ */
template <class R>
inline RetType & operator+=(const TensorStorage<T, ndim, R> & other) {
T * a = this->storage();
T * b = other.storage();
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
for (UInt i = 0; i < _size; ++i) {
*(a++) += *(b++);
}
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
template <class R>
inline RetType & operator-=(const TensorStorage<T, ndim, R> & other) {
T * a = this->storage();
T * b = other.storage();
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
for (UInt i = 0; i < _size; ++i) {
*(a++) -= *(b++);
}
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator+=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i) {
*(a++) += x;
}
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator-=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i) {
*(a++) -= x;
}
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
inline RetType & operator*=(const T & x) {
T * a = this->storage();
for (UInt i = 0; i < _size; ++i) {
*(a++) *= x;
}
return *(static_cast<RetType *>(this));
}
/* ---------------------------------------------------------------------- */
inline RetType & operator/=(const T & x) {
T * a = this->values;
for (UInt i = 0; i < _size; ++i) {
*(a++) /= x;
}
return *(static_cast<RetType *>(this));
}
/// \f[Y = \alpha X + Y\f]
inline RetType & aXplusY(const TensorStorage & other, const T & alpha = 1.) {
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be subtracted");
Math::aXplusY(this->_size, alpha, other.storage(), this->storage());
return *(static_cast<RetType *>(this));
}
/* ------------------------------------------------------------------------ */
T * storage() const { return values; }
UInt size() const { return _size; }
UInt size(UInt i) const {
AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
<< " dimensions, not "
<< (i + 1));
return n[i];
};
/* ------------------------------------------------------------------------ */
inline void set(const T & t) { std::fill_n(values, _size, t); };
inline void zero() { this->set(0.); };
-
+
template <class TensorType> inline void copy(const TensorType & other) {
AKANTU_DEBUG_ASSERT(
_size == other.size(),
"The two tensors do not have the same size, they cannot be copied");
std::copy_n(other.storage(), _size, values);
}
bool isWrapped() const { return this->wrapped; }
protected:
inline void computeSize() {
_size = 1;
for (UInt d = 0; d < ndim; ++d) {
_size *= this->n[d];
}
}
protected:
template <typename R, NormType norm_type> struct NormHelper {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it) {
_norm += std::pow(std::abs(*it), norm_type);
}
return std::pow(_norm, 1. / norm_type);
}
- };
+ }; // namespace akantu
template <typename R> struct NormHelper<R, L_1> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it) {
_norm += std::abs(*it);
}
return _norm;
}
};
template <typename R> struct NormHelper<R, L_2> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it) {
_norm += *it * *it;
}
return sqrt(_norm);
}
};
template <typename R> struct NormHelper<R, L_inf> {
template <class Ten> static R norm(const Ten & ten) {
R _norm = 0.;
R * it = ten.storage();
R * end = ten.storage() + ten.size();
for (; it < end; ++it) {
_norm = std::max(std::abs(*it), _norm);
}
return _norm;
}
};
public:
/*----------------------------------------------------------------------- */
/// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
/// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
/// @f]
template <NormType norm_type> inline T norm() const {
return NormHelper<T, norm_type>::norm(*this);
}
protected:
std::array<UInt, ndim> n{};
UInt _size{0};
T * values{nullptr};
bool wrapped{false};
};
/* -------------------------------------------------------------------------- */
/* Vector */
/* -------------------------------------------------------------------------- */
template <typename T> class Vector : public TensorStorage<T, 1, Vector<T>> {
using parent = TensorStorage<T, 1, Vector<T>>;
public:
using value_type = typename parent::value_type;
using proxy = VectorProxy<T>;
public:
Vector() : parent() {}
explicit Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) {}
Vector(T * data, UInt n) : parent(data, n, 0, 0) {}
Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) {}
Vector(const TensorProxy<T, 1, Vector> & src) : parent(src) {}
Vector(std::initializer_list<T> list) : parent(list.size(), 0, 0, T()) {
UInt i = 0;
for (auto val : list) {
operator()(i++) = val;
}
}
public:
using iterator = types::details::vector_iterator<T &>;
using const_iterator = types::details::vector_iterator<const T &>;
iterator begin() { return iterator(this->storage()); }
iterator end() { return iterator(this->storage() + this->size()); }
const_iterator begin() const { return const_iterator(this->storage()); }
const_iterator end() const {
return const_iterator(this->storage() + this->size());
}
public:
~Vector() override = default;
/* ------------------------------------------------------------------------ */
inline Vector & operator=(const Vector & src) {
parent::operator=(src);
return *this;
}
inline Vector & operator=(Vector && src) noexcept = default;
/* ------------------------------------------------------------------------ */
inline T & operator()(UInt i) {
AKANTU_DEBUG_ASSERT((i < this->n[0]),
"Access out of the vector! "
<< "Index (" << i
<< ") is out of the vector of size (" << this->n[0]
<< ")");
return *(this->values + i);
}
inline const T & operator()(UInt i) const {
AKANTU_DEBUG_ASSERT((i < this->n[0]),
"Access out of the vector! "
<< "Index (" << i
<< ") is out of the vector of size (" << this->n[0]
<< ")");
return *(this->values + i);
}
inline T & operator[](UInt i) { return this->operator()(i); }
inline const T & operator[](UInt i) const { return this->operator()(i); }
/* ------------------------------------------------------------------------ */
inline Vector<T> & operator*=(Real x) { return parent::operator*=(x); }
inline Vector<T> & operator/=(Real x) { return parent::operator/=(x); }
/* ------------------------------------------------------------------------ */
inline Vector<T> & operator*=(const Vector<T> & vect) {
AKANTU_DEBUG_ASSERT(this->_size == vect._size,
"The vectors have non matching sizes");
T * a = this->storage();
T * b = vect.storage();
for (UInt i = 0; i < this->_size; ++i) {
*(a++) *= *(b++);
}
return *this;
}
/* ------------------------------------------------------------------------ */
inline Real dot(const Vector<T> & vect) const {
return Math::vectorDot(this->values, vect.storage(), this->_size);
}
/* ------------------------------------------------------------------------ */
inline Real mean() const {
Real mean = 0;
T * a = this->storage();
for (UInt i = 0; i < this->_size; ++i) {
mean += *(a++);
}
return mean / this->_size;
}
/* ------------------------------------------------------------------------ */
inline Vector & crossProduct(const Vector<T> & v1, const Vector<T> & v2) {
AKANTU_DEBUG_ASSERT(this->size() == 3,
"crossProduct is only defined in 3D (n=" << this->size()
<< ")");
AKANTU_DEBUG_ASSERT(
this->size() == v1.size() && this->size() == v2.size(),
"crossProduct is not a valid operation non matching size vectors");
Math::vectorProduct3(v1.storage(), v2.storage(), this->values);
return *this;
}
inline Vector crossProduct(const Vector<T> & v) {
Vector<T> tmp(this->size());
tmp.crossProduct(*this, v);
return tmp;
}
/* ------------------------------------------------------------------------ */
inline void solve(const Matrix<T> & A, const Vector<T> & b) {
AKANTU_DEBUG_ASSERT(
this->size() == A.rows() && this->_size == A.cols(),
"The size of the solution vector mismatches the size of the matrix");
AKANTU_DEBUG_ASSERT(
this->_size == b._size,
"The rhs vector has a mismatch in size with the matrix");
Math::solve(this->_size, A.storage(), this->values, b.storage());
}
/* ------------------------------------------------------------------------ */
template <bool tr_A>
inline void mul(const Matrix<T> & A, const Vector<T> & x, T alpha = T(1));
/* ------------------------------------------------------------------------ */
inline Real norm() const { return parent::template norm<L_2>(); }
template <NormType nt> inline Real norm() const {
return parent::template norm<nt>();
}
/* ------------------------------------------------------------------------ */
inline Vector<Real> & normalize() {
Real n = norm();
operator/=(n);
return *this;
}
/* ------------------------------------------------------------------------ */
/// norm of (*this - x)
inline Real distance(const Vector<T> & y) const {
Real * vx = this->values;
Real * vy = y.storage();
Real sum_2 = 0;
for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy) { // NOLINT
sum_2 += (*vx - *vy) * (*vx - *vy);
}
return sqrt(sum_2);
}
/* ------------------------------------------------------------------------ */
inline bool equal(const Vector<T> & v,
Real tolerance = Math::getTolerance()) const {
T * a = this->storage();
T * b = v.storage();
UInt i = 0;
while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) {
++i;
}
return i == this->_size;
}
/* ------------------------------------------------------------------------ */
inline short compare(const Vector<T> & v,
Real tolerance = Math::getTolerance()) const {
T * a = this->storage();
T * b = v.storage();
for (UInt i(0); i < this->_size; ++i, ++a, ++b) {
if (std::abs(*a - *b) > tolerance) {
return (((*a - *b) > tolerance) ? 1 : -1);
}
}
return 0;
}
/* ------------------------------------------------------------------------ */
inline bool operator==(const Vector<T> & v) const { return equal(v); }
inline bool operator!=(const Vector<T> & v) const { return !operator==(v); }
inline bool operator<(const Vector<T> & v) const { return compare(v) == -1; }
inline bool operator>(const Vector<T> & v) const { return compare(v) == 1; }
template <typename Func, typename Acc>
decltype(auto) accumulate(const Vector<T> & v, Acc && accumulator,
Func && func) const {
T * a = this->storage();
T * b = v.storage();
for (UInt i(0); i < this->_size; ++i, ++a, ++b) {
accumulator = func(*a, *b, std::forward<Acc>(accumulator));
}
return accumulator;
}
inline bool operator<=(const Vector<T> & v) const {
bool res = true;
return accumulate(v, res, [](auto && a, auto && b, auto && accumulator) {
return accumulator & (a <= b);
});
}
inline bool operator>=(const Vector<T> & v) const {
bool res = true;
return accumulate(v, res, [](auto && a, auto && b, auto && accumulator) {
return accumulator & (a >= b);
});
}
/* ------------------------------------------------------------------------ */
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) {
;
}
stream << "[";
for (UInt i = 0; i < this->_size; ++i) {
if (i != 0) {
stream << ", ";
}
stream << this->values[i];
}
stream << "]";
}
/* ---------------------------------------------------------------------- */
static inline Vector<T> zeros(UInt n) {
Vector<T> tmp(n);
tmp.set(T());
return tmp;
}
};
using RVector = Vector<Real>;
/* ------------------------------------------------------------------------ */
template <>
inline bool Vector<UInt>::equal(const Vector<UInt> & v,
__attribute__((unused)) Real tolerance) const {
UInt * a = this->storage();
UInt * b = v.storage();
UInt i = 0;
while (i < this->_size && (*(a++) == *(b++))) {
++i;
}
return i == this->_size;
}
/* -------------------------------------------------------------------------- */
namespace types {
namespace details {
template <typename Mat> class column_iterator {
public:
using difference_type = std::ptrdiff_t;
using value_type = decltype(std::declval<Mat>().operator()(0));
using pointer = value_type *;
using reference = value_type &;
using iterator_category = std::input_iterator_tag;
column_iterator(Mat & mat, UInt col) : mat(mat), col(col) {}
decltype(auto) operator*() { return mat(col); }
decltype(auto) operator++() {
++col;
AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound");
return *this;
}
decltype(auto) operator++(int) {
auto tmp = *this;
++col;
AKANTU_DEBUG_ASSERT(col <= mat.cols(), "The iterator is out of bound");
return tmp;
}
bool operator!=(const column_iterator & other) const {
return col != other.col;
}
bool operator==(const column_iterator & other) const {
return not operator!=(other);
}
private:
Mat & mat;
UInt col;
};
} // namespace details
} // namespace types
/* ------------------------------------------------------------------------ */
/* Matrix */
/* ------------------------------------------------------------------------ */
template <typename T> class Matrix : public TensorStorage<T, 2, Matrix<T>> {
using parent = TensorStorage<T, 2, Matrix<T>>;
public:
using value_type = typename parent::value_type;
using proxy = MatrixProxy<T>;
public:
Matrix() : parent() {}
Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) {}
Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) {}
Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) {}
Matrix(const MatrixProxy<T> & src) : parent(src) {}
Matrix(std::initializer_list<std::initializer_list<T>> list) {
static_assert(std::is_trivially_copyable<T>{},
"Cannot create a tensor on non trivial types");
std::size_t n = 0;
std::size_t m = list.size();
for (auto row : list) {
n = std::max(n, row.size());
}
DimHelper<2>::setDims(m, n, 0, this->n);
this->computeSize();
this->values = new T[this->_size];
this->set(0);
UInt i{0};
UInt j{0};
for (auto & row : list) {
for (auto & val : row) {
at(i, j++) = val;
}
++i;
j = 0;
}
}
~Matrix() override = default;
/* ------------------------------------------------------------------------ */
inline Matrix & operator=(const Matrix & src) {
parent::operator=(src);
return *this;
}
inline Matrix & operator=(Matrix && src) noexcept = default;
+
public:
/* ---------------------------------------------------------------------- */
UInt rows() const { return this->n[0]; }
UInt cols() const { return this->n[1]; }
/* ---------------------------------------------------------------------- */
inline T & at(UInt i, UInt j) {
AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
"Access out of the matrix! "
<< "Index (" << i << ", " << j
<< ") is out of the matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return *(this->values + i + j * this->n[0]);
}
inline const T & at(UInt i, UInt j) const {
AKANTU_DEBUG_ASSERT(((i < this->n[0]) && (j < this->n[1])),
"Access out of the matrix! "
<< "Index (" << i << ", " << j
<< ") is out of the matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return *(this->values + i + j * this->n[0]);
}
/* ------------------------------------------------------------------------ */
inline T & operator()(UInt i, UInt j) { return this->at(i, j); }
inline const T & operator()(UInt i, UInt j) const { return this->at(i, j); }
/// give a line vector wrapped on the column i
inline VectorProxy<T> operator()(UInt j) {
AKANTU_DEBUG_ASSERT(j < this->n[1],
"Access out of the matrix! "
<< "You are trying to access the column vector "
<< j << " in a matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
}
inline VectorProxy<T> operator()(UInt j) const {
AKANTU_DEBUG_ASSERT(j < this->n[1],
"Access out of the matrix! "
<< "You are trying to access the column vector "
<< j << " in a matrix of size (" << this->n[0]
<< ", " << this->n[1] << ")");
return VectorProxy<T>(this->values + j * this->n[0], this->n[0]);
}
public:
decltype(auto) begin() {
return types::details::column_iterator<Matrix<T>>(*this, 0);
}
decltype(auto) begin() const {
return types::details::column_iterator<const Matrix<T>>(*this, 0);
}
decltype(auto) end() {
return types::details::column_iterator<Matrix<T>>(*this, this->cols());
}
decltype(auto) end() const {
return types::details::column_iterator<const Matrix<T>>(*this,
this->cols());
}
/* ------------------------------------------------------------------------ */
inline void block(const Matrix & block, UInt pos_i, UInt pos_j) {
AKANTU_DEBUG_ASSERT(pos_i + block.rows() <= rows(),
"The block size or position are not correct");
AKANTU_DEBUG_ASSERT(pos_i + block.cols() <= cols(),
"The block size or position are not correct");
for (UInt i = 0; i < block.rows(); ++i) {
for (UInt j = 0; j < block.cols(); ++j) {
this->at(i + pos_i, j + pos_j) = block(i, j);
}
}
}
inline Matrix block(UInt pos_i, UInt pos_j, UInt block_rows,
UInt block_cols) const {
AKANTU_DEBUG_ASSERT(pos_i + block_rows <= rows(),
"The block size or position are not correct");
AKANTU_DEBUG_ASSERT(pos_i + block_cols <= cols(),
"The block size or position are not correct");
Matrix block(block_rows, block_cols);
for (UInt i = 0; i < block_rows; ++i) {
for (UInt j = 0; j < block_cols; ++j) {
block(i, j) = this->at(i + pos_i, j + pos_j);
}
}
return block;
}
inline T & operator[](UInt idx) { return *(this->values + idx); };
inline const T & operator[](UInt idx) const { return *(this->values + idx); };
/* ---------------------------------------------------------------------- */
inline Matrix operator*(const Matrix & B) const {
Matrix C(this->rows(), B.cols());
C.mul<false, false>(*this, B);
return C;
}
/* ----------------------------------------------------------------------- */
inline Matrix & operator*=(const T & x) { return parent::operator*=(x); }
inline Matrix & operator*=(const Matrix & B) {
Matrix C(*this);
this->mul<false, false>(C, B);
return *this;
}
/* ---------------------------------------------------------------------- */
template <bool tr_A, bool tr_B>
inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) {
UInt k = A.cols();
if (tr_A) {
k = A.rows();
}
#ifndef AKANTU_NDEBUG
if (tr_B) {
AKANTU_DEBUG_ASSERT(k == B.cols(),
"matrices to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->cols() == B.rows(),
"matrices to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(k == B.rows(),
"matrices to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->cols() == B.cols(),
"matrices to multiply have no fit dimensions");
}
if (tr_A) {
AKANTU_DEBUG_ASSERT(this->rows() == A.cols(),
"matrices to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(this->rows() == A.rows(),
"matrices to multiply have no fit dimensions");
}
#endif // AKANTU_NDEBUG
Math::matMul<tr_A, tr_B>(this->rows(), this->cols(), k, alpha, A.storage(),
B.storage(), 0., this->storage());
}
/* ---------------------------------------------------------------------- */
inline void outerProduct(const Vector<T> & A, const Vector<T> & B) {
AKANTU_DEBUG_ASSERT(
A.size() == this->rows() && B.size() == this->cols(),
"A and B are not compatible with the size of the matrix");
for (UInt i = 0; i < this->rows(); ++i) {
for (UInt j = 0; j < this->cols(); ++j) {
this->values[i + j * this->rows()] += A[i] * B[j];
}
}
}
private:
class EigenSorter {
public:
EigenSorter(const Vector<T> & eigs) : eigs(eigs) {}
bool operator()(const UInt & a, const UInt & b) const {
return (eigs(a) > eigs(b));
}
private:
const Vector<T> & eigs;
};
public:
/* ---------------------------------------------------------------------- */
inline void eig(Vector<T> & eigenvalues, Matrix<T> & eigenvectors,
bool sort = true) const {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"eig is not a valid operation on a rectangular matrix");
AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(),
"eigenvalues should be of size " << this->cols()
<< ".");
#ifndef AKANTU_NDEBUG
if (eigenvectors.storage() != nullptr) {
AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) &&
(eigenvectors.rows() == this->cols()),
"Eigenvectors needs to be a square matrix of size "
<< this->cols() << " x " << this->cols() << ".");
}
#endif
Matrix<T> tmp = *this;
Vector<T> tmp_eigs(eigenvalues.size());
Matrix<T> tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols());
if (tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0) {
Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage());
} else {
Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(),
tmp_eig_vects.storage());
}
if (not sort) {
eigenvalues = tmp_eigs;
eigenvectors = tmp_eig_vects;
return;
}
Vector<UInt> perm(eigenvalues.size());
for (UInt i = 0; i < perm.size(); ++i) {
perm(i) = i;
}
std::sort(perm.storage(), perm.storage() + perm.size(),
EigenSorter(tmp_eigs));
for (UInt i = 0; i < perm.size(); ++i) {
eigenvalues(i) = tmp_eigs(perm(i));
}
if (tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0) {
for (UInt i = 0; i < perm.size(); ++i) {
for (UInt j = 0; j < eigenvectors.rows(); ++j) {
eigenvectors(j, i) = tmp_eig_vects(j, perm(i));
}
}
}
}
/* ---------------------------------------------------------------------- */
inline void eig(Vector<T> & eigenvalues) const {
Matrix<T> empty;
eig(eigenvalues, empty);
}
/* ---------------------------------------------------------------------- */
inline void eye(T alpha = 1.) {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"eye is not a valid operation on a rectangular matrix");
this->zero();
for (UInt i = 0; i < this->cols(); ++i) {
this->values[i + i * this->rows()] = alpha;
}
}
/* ---------------------------------------------------------------------- */
static inline Matrix<T> eye(UInt m, T alpha = 1.) {
Matrix<T> tmp(m, m);
tmp.eye(alpha);
return tmp;
}
/* ---------------------------------------------------------------------- */
inline T trace() const {
AKANTU_DEBUG_ASSERT(
this->cols() == this->rows(),
"trace is not a valid operation on a rectangular matrix");
T trace = 0.;
for (UInt i = 0; i < this->rows(); ++i) {
trace += this->values[i + i * this->rows()];
}
return trace;
}
/* ---------------------------------------------------------------------- */
inline Matrix transpose() const {
Matrix tmp(this->cols(), this->rows());
for (UInt i = 0; i < this->rows(); ++i) {
for (UInt j = 0; j < this->cols(); ++j) {
tmp(j, i) = operator()(i, j);
}
}
return tmp;
}
/* ---------------------------------------------------------------------- */
inline void inverse(const Matrix & A) {
AKANTU_DEBUG_ASSERT(A.cols() == A.rows(),
"inv is not a valid operation on a rectangular matrix");
AKANTU_DEBUG_ASSERT(this->cols() == A.cols(),
"the matrix should have the same size as its inverse");
if (this->cols() == 1) {
*this->values = 1. / *A.storage();
} else if (this->cols() == 2) {
Math::inv2(A.storage(), this->values);
} else if (this->cols() == 3) {
Math::inv3(A.storage(), this->values);
} else {
Math::inv(this->cols(), A.storage(), this->values);
}
}
inline Matrix inverse() {
Matrix inv(this->rows(), this->cols());
inv.inverse(*this);
return inv;
}
/* --------------------------------------------------------------------- */
inline T det() const {
AKANTU_DEBUG_ASSERT(this->cols() == this->rows(),
"inv is not a valid operation on a rectangular matrix");
if (this->cols() == 1) {
return *(this->values);
}
if (this->cols() == 2) {
return Math::det2(this->values);
}
if (this->cols() == 3) {
return Math::det3(this->values);
}
return Math::det(this->cols(), this->values);
}
/* --------------------------------------------------------------------- */
inline T doubleDot(const Matrix<T> & other) const {
AKANTU_DEBUG_ASSERT(
this->cols() == this->rows(),
"doubleDot is not a valid operation on a rectangular matrix");
if (this->cols() == 1) {
return *(this->values) * *(other.storage());
}
if (this->cols() == 2) {
return Math::matrixDoubleDot22(this->values, other.storage());
}
if (this->cols() == 3) {
return Math::matrixDoubleDot33(this->values, other.storage());
}
AKANTU_ERROR("doubleDot is not defined for other spatial dimensions"
<< " than 1, 2 or 3.");
}
/* ---------------------------------------------------------------------- */
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) {
;
}
stream << "[";
for (UInt i = 0; i < this->n[0]; ++i) {
if (i != 0) {
stream << ", ";
}
stream << "[";
for (UInt j = 0; j < this->n[1]; ++j) {
if (j != 0) {
stream << ", ";
}
stream << operator()(i, j);
}
stream << "]";
}
stream << "]";
};
};
/* ------------------------------------------------------------------------ */
template <typename T>
template <bool tr_A>
inline void Vector<T>::mul(const Matrix<T> & A, const Vector<T> & x, T alpha) {
#ifndef AKANTU_NDEBUG
UInt n = x.size();
if (tr_A) {
AKANTU_DEBUG_ASSERT(n == A.rows(),
"matrix and vector to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->size() == A.cols(),
"matrix and vector to multiply have no fit dimensions");
} else {
AKANTU_DEBUG_ASSERT(n == A.cols(),
"matrix and vector to multiply have no fit dimensions");
AKANTU_DEBUG_ASSERT(this->size() == A.rows(),
"matrix and vector to multiply have no fit dimensions");
}
#endif
Math::matVectMul<tr_A>(A.rows(), A.cols(), alpha, A.storage(), x.storage(),
0., this->storage());
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Matrix<T> & _this) {
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline std::ostream & operator<<(std::ostream & stream,
const Vector<T> & _this) {
_this.printself(stream);
return stream;
}
/* ------------------------------------------------------------------------ */
/* Tensor3 */
/* ------------------------------------------------------------------------ */
template <typename T> class Tensor3 : public TensorStorage<T, 3, Tensor3<T>> {
using parent = TensorStorage<T, 3, Tensor3<T>>;
public:
using value_type = typename parent::value_type;
using proxy = Tensor3Proxy<T>;
public:
Tensor3() : parent(){};
Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) {}
Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) {}
Tensor3(const Tensor3 & src, bool deep_copy = true)
: parent(src, deep_copy) {}
Tensor3(const proxy & src) : parent(src) {}
public:
/* ------------------------------------------------------------------------ */
inline Tensor3 & operator=(const Tensor3 & src) {
parent::operator=(src);
return *this;
}
/* ---------------------------------------------------------------------- */
inline T & operator()(UInt i, UInt j, UInt k) {
AKANTU_DEBUG_ASSERT(
(i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the element "
<< "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
<< this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
}
inline const T & operator()(UInt i, UInt j, UInt k) const {
AKANTU_DEBUG_ASSERT(
(i < this->n[0]) && (j < this->n[1]) && (k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the element "
<< "(" << i << ", " << j << ", " << k << ") in a tensor of size ("
<< this->n[0] << ", " << this->n[1] << ", " << this->n[2] << ")");
return *(this->values + (k * this->n[0] + i) * this->n[1] + j);
}
inline MatrixProxy<T> operator()(UInt k) {
AKANTU_DEBUG_ASSERT((k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the slice " << k
<< " in a tensor3 of size (" << this->n[0] << ", "
<< this->n[1] << ", " << this->n[2] << ")");
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline MatrixProxy<T> operator()(UInt k) const {
AKANTU_DEBUG_ASSERT((k < this->n[2]),
"Access out of the tensor3! "
<< "You are trying to access the slice " << k
<< " in a tensor3 of size (" << this->n[0] << ", "
<< this->n[1] << ", " << this->n[2] << ")");
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline MatrixProxy<T> operator[](UInt k) {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
inline MatrixProxy<T> operator[](UInt k) const {
return MatrixProxy<T>(this->values + k * this->n[0] * this->n[1],
this->n[0], this->n[1]);
}
};
/* -------------------------------------------------------------------------- */
// support operations for the creation of other vectors
/* -------------------------------------------------------------------------- */
template <typename T>
Vector<T> operator*(const T & scalar, const Vector<T> & a) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator*(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Vector<T> operator/(const Vector<T> & a, const T & scalar) {
Vector<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Vector<T> operator*(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r *= b;
return r;
}
template <typename T>
Vector<T> operator+(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r += b;
return r;
}
template <typename T>
Vector<T> operator-(const Vector<T> & a, const Vector<T> & b) {
Vector<T> r(a);
r -= b;
return r;
}
template <typename T>
Vector<T> operator*(const Matrix<T> & A, const Vector<T> & b) {
Vector<T> r(b.size());
r.template mul<false>(A, b);
return r;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Matrix<T> operator*(const T & scalar, const Matrix<T> & a) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator*(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r *= scalar;
return r;
}
template <typename T>
Matrix<T> operator/(const Matrix<T> & a, const T & scalar) {
Matrix<T> r(a);
r /= scalar;
return r;
}
template <typename T>
Matrix<T> operator+(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r += b;
return r;
}
template <typename T>
Matrix<T> operator-(const Matrix<T> & a, const Matrix<T> & b) {
Matrix<T> r(a);
r -= b;
return r;
}
} // namespace akantu
#include <iterator>
namespace std {
template <typename R>
struct iterator_traits<::akantu::types::details::vector_iterator<R>> {
protected:
using iterator = ::akantu::types::details::vector_iterator<R>;
public:
using iterator_category = typename iterator::iterator_category;
using value_type = typename iterator::value_type;
using difference_type = typename iterator::difference_type;
using pointer = typename iterator::pointer;
using reference = typename iterator::reference;
};
template <typename Mat>
struct iterator_traits<::akantu::types::details::column_iterator<Mat>> {
protected:
using iterator = ::akantu::types::details::column_iterator<Mat>;
public:
using iterator_category = typename iterator::iterator_category;
using value_type = typename iterator::value_type;
using difference_type = typename iterator::difference_type;
using pointer = typename iterator::pointer;
using reference = typename iterator::reference;
};
} // namespace std
#endif /* AKANTU_AKA_TYPES_HH_ */
diff --git a/src/io/mesh_io/mesh_io_msh.cc b/src/io/mesh_io/mesh_io_msh.cc
index 9549a66a8..c6c76f2d3 100644
--- a/src/io/mesh_io/mesh_io_msh.cc
+++ b/src/io/mesh_io/mesh_io_msh.cc
@@ -1,1123 +1,1123 @@
/**
* @file mesh_io_msh.cc
*
* @author Dana Christen <dana.christen@gmail.com>
* @author Mauro Corrado <mauro.corrado@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Tue Feb 20 2018
*
* @brief Read/Write for MSH files generated by gmsh
*
*
* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu 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.
*
* Akantu 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 Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -----------------------------------------------------------------------------
Version (Legacy) 1.0
$NOD
number-of-nodes
node-number x-coord y-coord z-coord
...
$ENDNOD
$ELM
number-of-elements
elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
...
$ENDELM
-----------------------------------------------------------------------------
Version 2.1
$MeshFormat
version-number file-type data-size
$EndMeshFormat
$Nodes
number-of-nodes
node-number x-coord y-coord z-coord
...
$EndNodes
$Elements
number-of-elements
elm-number elm-type number-of-tags < tag > ... node-number-list
...
$EndElements
$PhysicalNames
number-of-names
physical-dimension physical-number "physical-name"
...
$EndPhysicalNames
$NodeData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
node-number value ...
...
$EndNodeData
$ElementData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
elm-number value ...
...
$EndElementData
$ElementNodeData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
elm-number number-of-nodes-per-element value ...
...
$ElementEndNodeData
-----------------------------------------------------------------------------
elem-type
1: 2-node line.
2: 3-node triangle.
3: 4-node quadrangle.
4: 4-node tetrahedron.
5: 8-node hexahedron.
6: 6-node prism.
7: 5-node pyramid.
8: 3-node second order line
9: 6-node second order triangle
10: 9-node second order quadrangle
11: 10-node second order tetrahedron
12: 27-node second order hexahedron
13: 18-node second order prism
14: 14-node second order pyramid
15: 1-node point.
16: 8-node second order quadrangle
17: 20-node second order hexahedron
18: 15-node second order prism
19: 13-node second order pyramid
20: 9-node third order incomplete triangle
21: 10-node third order triangle
22: 12-node fourth order incomplete triangle
23: 15-node fourth order triangle
24: 15-node fifth order incomplete triangle
25: 21-node fifth order complete triangle
26: 4-node third order edge
27: 5-node fourth order edge
28: 6-node fifth order edge
29: 20-node third order tetrahedron
30: 35-node fourth order tetrahedron
31: 56-node fifth order tetrahedron
-------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "element_group.hh"
#include "mesh_io.hh"
#include "mesh_utils.hh"
#include "node_group.hh"
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Methods Implentations */
/* -------------------------------------------------------------------------- */
MeshIOMSH::MeshIOMSH() {
canReadSurface = true;
canReadExtendedData = true;
_msh_nodes_per_elem[_msh_not_defined] = 0;
_msh_nodes_per_elem[_msh_segment_2] = 2;
_msh_nodes_per_elem[_msh_triangle_3] = 3;
_msh_nodes_per_elem[_msh_quadrangle_4] = 4;
_msh_nodes_per_elem[_msh_tetrahedron_4] = 4;
_msh_nodes_per_elem[_msh_hexahedron_8] = 8;
_msh_nodes_per_elem[_msh_prism_1] = 6;
_msh_nodes_per_elem[_msh_pyramid_1] = 1;
_msh_nodes_per_elem[_msh_segment_3] = 3;
_msh_nodes_per_elem[_msh_triangle_6] = 6;
_msh_nodes_per_elem[_msh_quadrangle_9] = 9;
_msh_nodes_per_elem[_msh_tetrahedron_10] = 10;
_msh_nodes_per_elem[_msh_hexahedron_27] = 27;
_msh_nodes_per_elem[_msh_hexahedron_20] = 20;
_msh_nodes_per_elem[_msh_prism_18] = 18;
_msh_nodes_per_elem[_msh_prism_15] = 15;
_msh_nodes_per_elem[_msh_pyramid_14] = 14;
_msh_nodes_per_elem[_msh_point] = 1;
_msh_nodes_per_elem[_msh_quadrangle_8] = 8;
_msh_to_akantu_element_types[_msh_not_defined] = _not_defined;
_msh_to_akantu_element_types[_msh_segment_2] = _segment_2;
_msh_to_akantu_element_types[_msh_triangle_3] = _triangle_3;
_msh_to_akantu_element_types[_msh_quadrangle_4] = _quadrangle_4;
_msh_to_akantu_element_types[_msh_tetrahedron_4] = _tetrahedron_4;
_msh_to_akantu_element_types[_msh_hexahedron_8] = _hexahedron_8;
_msh_to_akantu_element_types[_msh_prism_1] = _pentahedron_6;
_msh_to_akantu_element_types[_msh_pyramid_1] = _not_defined;
_msh_to_akantu_element_types[_msh_segment_3] = _segment_3;
_msh_to_akantu_element_types[_msh_triangle_6] = _triangle_6;
_msh_to_akantu_element_types[_msh_quadrangle_9] = _not_defined;
_msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10;
_msh_to_akantu_element_types[_msh_hexahedron_27] = _not_defined;
_msh_to_akantu_element_types[_msh_hexahedron_20] = _hexahedron_20;
_msh_to_akantu_element_types[_msh_prism_18] = _not_defined;
_msh_to_akantu_element_types[_msh_prism_15] = _pentahedron_15;
_msh_to_akantu_element_types[_msh_pyramid_14] = _not_defined;
_msh_to_akantu_element_types[_msh_point] = _point_1;
_msh_to_akantu_element_types[_msh_quadrangle_8] = _quadrangle_8;
_akantu_to_msh_element_types[_not_defined] = _msh_not_defined;
_akantu_to_msh_element_types[_segment_2] = _msh_segment_2;
_akantu_to_msh_element_types[_segment_3] = _msh_segment_3;
_akantu_to_msh_element_types[_triangle_3] = _msh_triangle_3;
_akantu_to_msh_element_types[_triangle_6] = _msh_triangle_6;
_akantu_to_msh_element_types[_tetrahedron_4] = _msh_tetrahedron_4;
_akantu_to_msh_element_types[_tetrahedron_10] = _msh_tetrahedron_10;
_akantu_to_msh_element_types[_quadrangle_4] = _msh_quadrangle_4;
_akantu_to_msh_element_types[_quadrangle_8] = _msh_quadrangle_8;
_akantu_to_msh_element_types[_hexahedron_8] = _msh_hexahedron_8;
_akantu_to_msh_element_types[_hexahedron_20] = _msh_hexahedron_20;
_akantu_to_msh_element_types[_pentahedron_6] = _msh_prism_1;
_akantu_to_msh_element_types[_pentahedron_15] = _msh_prism_15;
_akantu_to_msh_element_types[_point_1] = _msh_point;
#if defined(AKANTU_STRUCTURAL_MECHANICS)
_akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2;
_akantu_to_msh_element_types[_bernoulli_beam_3] = _msh_segment_2;
_akantu_to_msh_element_types[_discrete_kirchhoff_triangle_18] =
_msh_triangle_3;
#endif
std::map<ElementType, MSHElementType>::iterator it;
for (it = _akantu_to_msh_element_types.begin();
it != _akantu_to_msh_element_types.end(); ++it) {
UInt nb_nodes = _msh_nodes_per_elem[it->second];
std::vector<UInt> tmp(nb_nodes);
for (UInt i = 0; i < nb_nodes; ++i) {
tmp[i] = i;
}
switch (it->first) {
case _tetrahedron_10:
tmp[8] = 9;
tmp[9] = 8;
break;
case _pentahedron_6:
tmp[0] = 2;
tmp[1] = 0;
tmp[2] = 1;
tmp[3] = 5;
tmp[4] = 3;
tmp[5] = 4;
break;
case _pentahedron_15:
tmp[0] = 2;
tmp[1] = 0;
tmp[2] = 1;
tmp[3] = 5;
tmp[4] = 3;
tmp[5] = 4;
tmp[6] = 8;
tmp[8] = 11;
tmp[9] = 6;
tmp[10] = 9;
tmp[11] = 10;
tmp[12] = 14;
tmp[14] = 12;
break;
case _hexahedron_20:
tmp[9] = 11;
tmp[10] = 12;
tmp[11] = 9;
tmp[12] = 13;
tmp[13] = 10;
tmp[17] = 19;
tmp[18] = 17;
tmp[19] = 18;
break;
default:
// nothing to change
break;
}
_read_order[it->first] = tmp;
}
}
/* -------------------------------------------------------------------------- */
MeshIOMSH::~MeshIOMSH() = default;
/* -------------------------------------------------------------------------- */
namespace {
struct File {
std::string filename;
std::ifstream infile;
std::string line;
size_t current_line{0};
size_t first_node_number{std::numeric_limits<UInt>::max()},
last_node_number{0};
bool has_physical_names{false};
std::unordered_map<size_t, size_t> node_tags;
std::unordered_map<size_t, Element> element_tags;
double version{0};
int size_of_size_t{0};
Mesh & mesh;
MeshAccessor mesh_accessor;
std::multimap<std::pair<int, int>, int> entity_tag_to_physical_tags;
File(const std::string & filename, Mesh & mesh)
: filename(filename), mesh(mesh), mesh_accessor(mesh) {
infile.open(filename.c_str());
if (not infile.good()) {
AKANTU_EXCEPTION("Cannot open file " << filename);
}
}
~File() { infile.close(); }
auto good() { return infile.good(); }
std::stringstream get_line() {
std::string tmp_str;
if (infile.eof()) {
AKANTU_EXCEPTION("Reached the end of the file " << filename);
}
std::getline(infile, tmp_str);
line = trim(tmp_str);
++current_line;
return std::stringstream(line);
}
template <typename... Ts> void read_line(Ts &&... ts) {
auto && sstr = get_line();
(void)std::initializer_list<int>{
(sstr >> std::forward<decltype(ts)>(ts), 0)...};
}
};
} // namespace
/* -------------------------------------------------------------------------- */
template <typename File, typename Readers>
void MeshIOMSH::populateReaders2(File & file, Readers & readers) {
readers["$NOD"] = readers["$Nodes"] = [&](const std::string & /*unused*/) {
UInt nb_nodes;
file.read_line(nb_nodes);
Array<Real> & nodes = file.mesh_accessor.getNodes();
nodes.resize(nb_nodes);
file.mesh_accessor.setNbGlobalNodes(nb_nodes);
size_t index;
Vector<double> coord(3);
/// for each node, read the coordinates
for (auto && data : enumerate(make_view(nodes, nodes.getNbComponent()))) {
file.read_line(index, coord(0), coord(1), coord(2));
if (index > std::numeric_limits<UInt>::max()) {
AKANTU_EXCEPTION(
"There are more nodes in this files than the index type of akantu "
"can handle, consider recompiling with a bigger index type");
}
file.first_node_number = std::min(file.first_node_number, index);
file.last_node_number = std::max(file.last_node_number, index);
for (auto && coord_data : zip(std::get<1>(data), coord)) {
std::get<0>(coord_data) = std::get<1>(coord_data);
}
file.node_tags[index] = std::get<0>(data);
}
};
readers["$ELM"] = readers["$Elements"] = [&](const std::string & /*unused*/) {
UInt nb_elements;
file.read_line(nb_elements);
Int index;
UInt msh_type;
ElementType akantu_type;
for (UInt i = 0; i < nb_elements; ++i) {
auto && sstr_elem = file.get_line();
sstr_elem >> index;
sstr_elem >> msh_type;
/// get the connectivity vector depending on the element type
akantu_type =
this->_msh_to_akantu_element_types[MSHElementType(msh_type)];
if (akantu_type == _not_defined) {
AKANTU_DEBUG_WARNING("Unsuported element kind "
<< msh_type << " at line " << file.current_line);
continue;
}
Element elem{akantu_type, 0, _not_ghost};
auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type);
auto node_per_element = connectivity.getNbComponent();
auto & read_order = this->_read_order[akantu_type];
/// read tags informations
if (file.version < 2) {
Int tag0;
Int tag1;
Int nb_nodes; // reg-phys, reg-elem, number-of-nodes
sstr_elem >> tag0 >> tag1 >> nb_nodes;
auto & data0 =
file.mesh_accessor.template getData<UInt>("tag_0", akantu_type);
data0.push_back(tag0);
auto & data1 =
file.mesh_accessor.template getData<UInt>("tag_1", akantu_type);
data1.push_back(tag1);
} else if (file.version < 4) {
UInt nb_tags;
sstr_elem >> nb_tags;
for (UInt j = 0; j < nb_tags; ++j) {
Int tag;
sstr_elem >> tag;
auto & data = file.mesh_accessor.template getData<UInt>(
"tag_" + std::to_string(j), akantu_type);
data.push_back(tag);
}
}
Vector<UInt> local_connect(node_per_element);
for (UInt j = 0; j < node_per_element; ++j) {
UInt node_index;
sstr_elem >> node_index;
AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number,
"Node number not in range : line "
<< file.current_line);
local_connect(read_order[j]) = file.node_tags[node_index];
}
connectivity.push_back(local_connect);
elem.element = connectivity.size() - 1;
file.element_tags[index] = elem;
}
};
readers["$Periodic"] = [&](const std::string & /*unused*/) {
UInt nb_periodic_entities;
file.read_line(nb_periodic_entities);
file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(),
NodeFlag::_normal);
for (UInt p = 0; p < nb_periodic_entities; ++p) {
// dimension slave-tag master-tag
UInt dimension;
file.read_line(dimension);
// transformation
file.get_line();
// nb nodes
UInt nb_nodes;
file.read_line(nb_nodes);
for (UInt n = 0; n < nb_nodes; ++n) {
// slave master
auto && sstr = file.get_line();
// The info in the mesh seem inconsistent so they are ignored for now.
continue;
if (dimension == file.mesh.getSpatialDimension() - 1) {
UInt slave;
UInt master;
sstr >> slave;
sstr >> master;
file.mesh_accessor.addPeriodicSlave(file.node_tags[slave],
file.node_tags[master]);
}
}
}
// mesh_accessor.markMeshPeriodic();
};
}
/* -------------------------------------------------------------------------- */
template <typename File, typename Readers>
void MeshIOMSH::populateReaders4(File & file, Readers & readers) {
static std::map<int, std::string> entity_type{
{0, "points"},
{1, "curve"},
{2, "surface"},
{3, "volume"},
};
readers["$Entities"] = [&](const std::string & /*unused*/) {
size_t num_entity[4];
file.read_line(num_entity[0], num_entity[1], num_entity[2], num_entity[3]);
for (auto entity_dim : arange(4)) {
for (auto _ [[gnu::unused]] : arange(num_entity[entity_dim])) {
auto && sstr = file.get_line();
int tag;
double min_x;
double min_y;
double min_z;
double max_x;
double max_y;
double max_z;
size_t num_physical_tags;
sstr >> tag >> min_x >> min_y >> min_z;
if (entity_dim > 0 or file.version < 4.1) {
sstr >> max_x >> max_y >> max_z;
}
sstr >> num_physical_tags;
for (auto _ [[gnu::unused]] : arange(num_physical_tags)) {
int phys_tag;
sstr >> phys_tag;
std::string physical_name;
if (this->physical_names.find(phys_tag) ==
this->physical_names.end()) {
physical_name = "msh_block_" + std::to_string(phys_tag);
} else {
physical_name = this->physical_names[phys_tag];
}
if (not file.mesh.elementGroupExists(physical_name)) {
file.mesh.createElementGroup(physical_name, entity_dim);
} else {
file.mesh.getElementGroup(physical_name).addDimension(entity_dim);
}
file.entity_tag_to_physical_tags.insert(
std::make_pair(std::make_pair(tag, entity_dim), phys_tag));
}
}
}
};
readers["$Nodes"] = [&](const std::string & /*unused*/) {
size_t num_blocks;
size_t num_nodes;
if (file.version >= 4.1) {
file.read_line(num_blocks, num_nodes, file.first_node_number,
file.last_node_number);
} else {
file.read_line(num_blocks, num_nodes);
}
auto & nodes = file.mesh_accessor.getNodes();
nodes.reserve(num_nodes);
file.mesh_accessor.setNbGlobalNodes(num_nodes);
if (num_nodes > std::numeric_limits<UInt>::max()) {
AKANTU_EXCEPTION(
"There are more nodes in this files than the index type of akantu "
"can handle, consider recompiling with a bigger index type");
}
size_t node_id{0};
for (auto block [[gnu::unused]] : arange(num_blocks)) {
int entity_dim;
int entity_tag;
int parametric;
size_t num_nodes_in_block;
Vector<double> pos(3);
Vector<double> real_pos(nodes.getNbComponent());
if (file.version >= 4.1) {
file.read_line(entity_dim, entity_tag, parametric, num_nodes_in_block);
if (parametric) {
AKANTU_EXCEPTION(
"Akantu does not support parametric nodes in msh files");
}
for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
size_t tag;
file.read_line(tag);
file.node_tags[tag] = node_id;
++node_id;
}
for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
file.read_line(pos(_x), pos(_y), pos(_z));
for (auto && data : zip(real_pos, pos)) {
std::get<0>(data) = std::get<1>(data);
}
nodes.push_back(real_pos);
}
} else {
file.read_line(entity_tag, entity_dim, parametric, num_nodes_in_block);
for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
size_t tag;
file.read_line(tag, pos(_x), pos(_y), pos(_z));
if (file.version < 4.1) {
file.first_node_number = std::min(file.first_node_number, tag);
file.last_node_number = std::max(file.last_node_number, tag);
}
for (auto && data : zip(real_pos, pos)) {
std::get<0>(data) = std::get<1>(data);
}
nodes.push_back(real_pos);
file.node_tags[tag] = node_id;
++node_id;
}
}
}
};
readers["$Elements"] = [&](const std::string & /*unused*/) {
size_t num_blocks;
size_t num_elements;
file.read_line(num_blocks, num_elements);
for (auto block [[gnu::unused]] : arange(num_blocks)) {
int entity_dim;
int entity_tag;
int element_type;
size_t num_elements_in_block;
if (file.version >= 4.1) {
file.read_line(entity_dim, entity_tag, element_type,
num_elements_in_block);
} else {
file.read_line(entity_tag, entity_dim, element_type,
num_elements_in_block);
}
/// get the connectivity vector depending on the element type
auto && akantu_type =
this->_msh_to_akantu_element_types[(MSHElementType)element_type];
if (akantu_type == _not_defined) {
AKANTU_DEBUG_WARNING("Unsuported element kind " << element_type
<< " at line "
<< file.current_line);
continue;
}
Element elem{akantu_type, 0, _not_ghost};
auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type);
Vector<UInt> local_connect(connectivity.getNbComponent());
auto && read_order = this->_read_order[akantu_type];
auto & data0 =
file.mesh_accessor.template getData<UInt>("tag_0", akantu_type);
data0.resize(data0.size() + num_elements_in_block, 0);
auto & physical_data = file.mesh_accessor.template getData<std::string>(
"physical_names", akantu_type);
physical_data.resize(physical_data.size() + num_elements_in_block, "");
for (auto _ [[gnu::unused]] : arange(num_elements_in_block)) {
auto && sstr_elem = file.get_line();
size_t elem_tag;
sstr_elem >> elem_tag;
for (auto && c : arange(connectivity.getNbComponent())) {
size_t node_tag;
sstr_elem >> node_tag;
AKANTU_DEBUG_ASSERT(node_tag <= file.last_node_number,
"Node number not in range : line "
<< file.current_line);
node_tag = file.node_tags[node_tag];
local_connect(read_order[c]) = node_tag;
}
connectivity.push_back(local_connect);
elem.element = connectivity.size() - 1;
file.element_tags[elem_tag] = elem;
auto range = file.entity_tag_to_physical_tags.equal_range(
std::make_pair(entity_tag, entity_dim));
bool first = true;
for (auto it = range.first; it != range.second; ++it) {
auto phys_it = this->physical_names.find(it->second);
if (first) {
data0(elem.element) =
it->second; // for compatibility with version 2
if (phys_it != this->physical_names.end()) {
physical_data(elem.element) = phys_it->second;
}
first = false;
}
if (phys_it != this->physical_names.end()) {
file.mesh.getElementGroup(phys_it->second).add(elem, true, false);
}
}
}
}
for (auto && element_group : file.mesh.iterateElementGroups()) {
element_group.getNodeGroup().optimize();
}
};
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::read(const std::string & filename, Mesh & mesh) {
File file(filename, mesh);
std::map<std::string, std::function<void(const std::string &)>> readers;
readers["$MeshFormat"] = [&](const std::string & /*unused*/) {
auto && sstr = file.get_line();
int format;
sstr >> file.version >> format;
if (format != 0) {
AKANTU_ERROR("This reader can only read ASCII files.");
}
if (file.version > 2) {
sstr >> file.size_of_size_t;
if (file.size_of_size_t > int(sizeof(UInt))) {
AKANTU_DEBUG_INFO("The size of the indexes in akantu might be to small "
"to read this file (akantu "
<< sizeof(UInt) << " vs. msh file "
<< file.size_of_size_t << ")");
}
}
if (file.version < 4) {
this->populateReaders2(file, readers);
} else {
this->populateReaders4(file, readers);
}
};
auto && read_data = [&](auto && entity_tags, auto && get_data,
auto && read_data) {
auto read_data_tags = [&](auto x) {
UInt number_of_tags{0};
file.read_line(number_of_tags);
std::vector<decltype(x)> tags(number_of_tags);
for (auto && tag : tags) {
file.read_line(tag);
}
return tags;
};
auto && string_tags = read_data_tags(std::string{});
auto && real_tags [[gnu::unused]] = read_data_tags(double{});
auto && int_tags = read_data_tags(int{});
for (auto & s : string_tags) {
s = trim(s, '"');
}
auto id = string_tags[0];
auto size = int_tags[2];
auto nb_component = int_tags[1];
auto & data = get_data(id, size, nb_component);
for (auto n [[gnu::unused]] : arange(size)) {
auto && sstr = file.get_line();
size_t tag;
sstr >> tag;
const auto & entity = entity_tags[tag];
read_data(entity, sstr, data, nb_component);
}
};
readers["$NodeData"] = [&](const std::string & /*unused*/) {
/* $NodeData
numStringTags(ASCII int)
stringTag(string) ...
numRealTags(ASCII int)
realTag(ASCII double) ...
numIntegerTags(ASCII int)
integerTag(ASCII int) ...
nodeTag(size_t) value(double) ...
...
$EndNodeData */
read_data(
file.node_tags,
[&](auto && id, auto && size [[gnu::unused]],
auto && nb_component [[gnu::unused]]) -> Array<double> & {
auto & data =
file.mesh.template getNodalData<double>(id, nb_component);
data.resize(size);
return data;
},
[&](auto && node, auto && sstr, auto && data,
auto && nb_component [[gnu::unused]]) {
for (auto c : arange(nb_component)) {
sstr >> data(node, c);
}
});
};
readers["$ElementData"] = [&](const std::string & /*unused*/) {
/* $ElementData
numStringTags(ASCII int)
stringTag(string) ...
numRealTags(ASCII int)
realTag(ASCII double) ...
numIntegerTags(ASCII int)
integerTag(ASCII int) ...
elementTag(size_t) value(double) ...
...
$EndElementData
*/
read_data(
file.element_tags,
[&](auto && id, auto && size [[gnu::unused]],
auto && nb_component
[[gnu::unused]]) -> ElementTypeMapArray<double> & {
file.mesh.template getElementalData<double>(id);
return file.mesh.template getElementalData<double>(id);
},
[&](auto && element, auto && sstr, auto && data, auto && nb_component) {
if (not data.exists(element.type)) {
data.alloc(mesh.getNbElement(element.type), nb_component,
element.type, element.ghost_type);
}
auto & data_array = data(element.type);
for (auto c : arange(nb_component)) {
sstr >> data_array(element.element, c);
}
});
};
readers["$ElementNodeData"] = [&](const std::string & /*unused*/) {
/* $ElementNodeData
numStringTags(ASCII int)
stringTag(string) ...
numRealTags(ASCII int)
realTag(ASCII double) ...
numIntegerTags(ASCII int)
integerTag(ASCII int) ...
elementTag(size_t) value(double) ...
...
$EndElementNodeData
*/
read_data(
file.element_tags,
[&](auto && id, auto && size [[gnu::unused]],
auto && nb_component
[[gnu::unused]]) -> ElementTypeMapArray<double> & {
file.mesh.template getElementalData<double>(id);
auto & data = file.mesh.template getElementalData<double>(id);
data.isNodal(true);
return data;
},
[&](auto && element, auto && sstr, auto && data, auto && nb_component) {
int nb_nodes_per_element;
sstr >> nb_nodes_per_element;
if (not data.exists(element.type)) {
data.alloc(mesh.getNbElement(element.type),
nb_component * nb_nodes_per_element, element.type,
element.ghost_type);
}
auto & data_array = data(element.type);
for (auto c : arange(nb_component)) {
sstr >> data_array(element.element, c);
}
});
};
readers["$PhysicalNames"] = [&](const std::string & /*unused*/) {
file.has_physical_names = true;
int num_of_phys_names;
file.read_line(num_of_phys_names); /// the format line
for (auto k [[gnu::unused]] : arange(num_of_phys_names)) {
int phys_name_id;
int phys_dim;
std::string phys_name;
file.read_line(phys_dim, phys_name_id, std::quoted(phys_name));
this->physical_names[phys_name_id] = phys_name;
}
};
readers["Unsupported"] = [&](const std::string & _block) {
std::string block = _block.substr(1);
AKANTU_DEBUG_WARNING("Unsupported block_kind " << block << " at line "
<< file.current_line);
auto && end_block = "$End" + block;
while (file.line != end_block) {
file.get_line();
}
};
while (file.good()) {
std::string block;
file.read_line(block);
auto && it = readers.find(block);
if (it != readers.end()) {
it->second(block);
std::string end_block;
file.read_line(end_block);
block = block.substr(1);
if (end_block != "$End" + block) {
AKANTU_EXCEPTION("The reader failed to properly read the block "
<< block << ". Expected a $End" << block << " at line "
<< file.current_line);
}
} else if (not block.empty()) {
readers["Unsupported"](block);
}
}
// mesh.updateTypesOffsets(_not_ghost);
if (file.version < 4) {
this->constructPhysicalNames("tag_0", mesh);
if (file.has_physical_names) {
mesh.createGroupsFromMeshData<std::string>("physical_names");
}
}
MeshUtils::fillElementToSubElementsData(mesh);
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) {
std::ofstream outfile;
const Array<Real> & nodes = mesh.getNodes();
outfile.open(filename.c_str());
outfile << "$MeshFormat"
<< "\n";
outfile << "2.2 0 8"
<< "\n";
outfile << "$EndMeshFormat"
<< "\n";
outfile << std::setprecision(std::numeric_limits<Real>::digits10);
outfile << "$Nodes"
<< "\n";
outfile << nodes.size() << "\n";
outfile << std::uppercase;
for (UInt i = 0; i < nodes.size(); ++i) {
Int offset = i * nodes.getNbComponent();
outfile << i + 1;
for (UInt j = 0; j < nodes.getNbComponent(); ++j) {
outfile << " " << nodes.storage()[offset + j];
}
for (UInt p = nodes.getNbComponent(); p < 3; ++p) {
outfile << " " << 0.;
}
outfile << "\n";
;
}
outfile << std::nouppercase;
outfile << "$EndNodes"
<< "\n";
outfile << "$Elements"
<< "\n";
Int nb_elements = 0;
for (auto && type :
mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) {
const Array<UInt> & connectivity = mesh.getConnectivity(type, _not_ghost);
nb_elements += connectivity.size();
}
outfile << nb_elements << "\n";
std::map<Element, size_t> element_to_msh_element;
UInt element_idx = 1;
- Element element;
+ auto element = ElementNull;
for (auto && type :
mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) {
const auto & connectivity = mesh.getConnectivity(type, _not_ghost);
element.type = type;
UInt * tag[2] = {nullptr, nullptr};
if (mesh.hasData<UInt>("tag_0", type, _not_ghost)) {
const auto & data_tag_0 = mesh.getData<UInt>("tag_0", type, _not_ghost);
tag[0] = data_tag_0.storage();
}
if (mesh.hasData<UInt>("tag_1", type, _not_ghost)) {
const auto & data_tag_1 = mesh.getData<UInt>("tag_1", type, _not_ghost);
tag[1] = data_tag_1.storage();
}
for (auto && data :
enumerate(make_view(connectivity, connectivity.getNbComponent()))) {
element.element = std::get<0>(data);
const auto & conn = std::get<1>(data);
- element_to_msh_element[element] = element_idx;
+ element_to_msh_element.insert(std::make_pair(element, element_idx));
outfile << element_idx << " " << _akantu_to_msh_element_types[type]
<< " 2";
/// \todo write the real data in the file
for (UInt t = 0; t < 2; ++t) {
if (tag[t] != nullptr) {
outfile << " " << tag[t][element.element];
} else {
outfile << " 0";
}
}
for (auto && c : conn) {
outfile << " " << c + 1;
}
outfile << "\n";
element_idx++;
}
}
outfile << "$EndElements"
<< "\n";
if (mesh.hasData(MeshDataType::_nodal)) {
auto && tags = mesh.getTagNames();
for (auto && tag : tags) {
auto type = mesh.getTypeCode(tag, MeshDataType::_nodal);
if (type != MeshDataTypeCode::_real) {
AKANTU_DEBUG_WARNING(
"The field "
<< tag << " is ignored by the MSH writer, msh files do not support "
<< type << " data");
continue;
}
auto && data = mesh.getNodalData<double>(tag);
outfile << "$NodeData"
<< "\n";
outfile << "1"
<< "\n";
outfile << "\"" << tag << "\"\n";
outfile << "1\n0.0"
<< "\n";
outfile << "3\n0"
<< "\n";
outfile << data.getNbComponent() << "\n";
outfile << data.size() << "\n";
for (auto && d : enumerate(make_view(data, data.getNbComponent()))) {
outfile << std::get<0>(d) + 1;
for (auto && v : std::get<1>(d)) {
outfile << " " << v;
}
outfile << "\n";
}
outfile << "$EndNodeData"
<< "\n";
}
}
if (mesh.hasData(MeshDataType::_elemental)) {
auto && tags = mesh.getTagNames();
for (auto && tag : tags) {
auto && data = mesh.getElementalData<double>(tag);
auto type = mesh.getTypeCode(tag, MeshDataType::_elemental);
if (type != MeshDataTypeCode::_real) {
AKANTU_DEBUG_WARNING(
"The field "
<< tag << " is ignored by the MSH writer, msh files do not support "
<< type << " data");
continue;
}
if (data.isNodal()) {
continue;
}
auto size = data.size();
if (size == 0) {
continue;
}
auto && nb_components = data.getNbComponents();
auto nb_component = nb_components(*(data.elementTypes().begin()));
outfile << "$ElementData"
<< "\n";
outfile << "1"
<< "\n";
outfile << "\"" << tag << "\"\n";
outfile << "1\n0.0"
<< "\n";
outfile << "3\n0"
<< "\n";
outfile << nb_component << "\n";
outfile << size << "\n";
Element element;
for (auto type : data.elementTypes()) {
element.type = type;
for (auto && _ :
enumerate(make_view(data(type), nb_components(type)))) {
element.element = std::get<0>(_);
outfile << element_to_msh_element[element];
for (auto && v : std::get<1>(_)) {
outfile << " " << v;
}
outfile << "\n";
}
}
outfile << "$EndElementData"
<< "\n";
}
}
outfile.close();
}
/* --------------------------------------------------------------------------
*/
} // namespace akantu