diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake
index c64a05a82..27a98df70 100644
--- a/packages/solid_mechanics.cmake
+++ b/packages/solid_mechanics.cmake
@@ -1,128 +1,129 @@
#===============================================================================
# @file solid_mechanics.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Mon Nov 21 2011
# @date last modification: Mon Jan 18 2016
#
# @brief package description for core
#
# @section LICENSE
#
# Copyright (©) 2010-2012, 2014, 2015 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/>.
#
#===============================================================================
package_declare(solid_mechanics DEFAULT ON
DESCRIPTION "Solid mechanics model"
DEPENDS core
)
package_declare_sources(solid_mechanics
model/solid_mechanics/material.cc
model/solid_mechanics/material.hh
model/solid_mechanics/material_inline_impl.cc
model/solid_mechanics/material_selector.hh
model/solid_mechanics/material_selector_tmpl.hh
model/solid_mechanics/materials/internal_field.hh
model/solid_mechanics/materials/internal_field_tmpl.hh
model/solid_mechanics/materials/random_internal_field.hh
model/solid_mechanics/materials/random_internal_field_tmpl.hh
model/solid_mechanics/solid_mechanics_model.cc
model/solid_mechanics/solid_mechanics_model.hh
model/solid_mechanics/solid_mechanics_model_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_io.cc
model/solid_mechanics/solid_mechanics_model_mass.cc
model/solid_mechanics/solid_mechanics_model_material.cc
model/solid_mechanics/solid_mechanics_model_tmpl.hh
model/solid_mechanics/solid_mechanics_model_event_handler.hh
model/solid_mechanics/materials/plane_stress_toolbox.hh
model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh
model/solid_mechanics/materials/material_core_includes.hh
model/solid_mechanics/materials/material_elastic.cc
model/solid_mechanics/materials/material_elastic.hh
model/solid_mechanics/materials/material_elastic_inline_impl.cc
model/solid_mechanics/materials/material_thermal.cc
model/solid_mechanics/materials/material_thermal.hh
model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
+ model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.cc
model/solid_mechanics/materials/material_elastic_orthotropic.cc
model/solid_mechanics/materials/material_elastic_orthotropic.hh
model/solid_mechanics/materials/material_damage/material_damage.hh
model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh
model/solid_mechanics/materials/material_damage/material_marigo.cc
model/solid_mechanics/materials/material_damage/material_marigo.hh
model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc
model/solid_mechanics/materials/material_damage/material_mazars.cc
model/solid_mechanics/materials/material_damage/material_mazars.hh
model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc
model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc
model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh
model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc
model/solid_mechanics/materials/material_plastic/material_plastic.cc
model/solid_mechanics/materials/material_plastic/material_plastic.hh
model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc
model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc
model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh
model/solid_mechanics/materials/material_non_local.hh
model/solid_mechanics/materials/material_non_local_tmpl.hh
model/solid_mechanics/materials/material_non_local_includes.hh
)
package_declare_material_infos(solid_mechanics
LIST AKANTU_CORE_MATERIAL_LIST
INCLUDE material_core_includes.hh
)
package_declare_documentation_files(solid_mechanics
manual-solidmechanicsmodel.tex
manual-constitutive-laws.tex
manual-lumping.tex
manual-appendix-materials.tex
figures/dynamic_analysis.png
figures/explicit_dynamic.pdf
figures/explicit_dynamic.svg
figures/static.pdf
figures/static.svg
figures/hooke_law.pdf
figures/implicit_dynamic.pdf
figures/implicit_dynamic.svg
figures/problemDomain.pdf_tex
figures/problemDomain.pdf
figures/static_analysis.png
figures/stress_strain_el.pdf
figures/tangent.pdf
figures/tangent.svg
figures/stress_strain_neo.pdf
figures/visco_elastic_law.pdf
figures/isotropic_hardening_plasticity.pdf
figures/stress_strain_visco.pdf
)
package_declare_extra_files_to_package(solid_mechanics
SOURCES
model/solid_mechanics/material_list.hh.in
)
diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index 62d998fde..203d4d318 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,1302 +1,1302 @@
/**
* @file aka_types.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Feb 17 2011
* @date last modification: Fri Jan 22 2016
*
* @brief description of the "simple" types
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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_error.hh"
#include "aka_fwd.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, UInt dims[dim]);
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<1> {
static inline void setDims(UInt m, __attribute__((unused)) UInt n,
__attribute__((unused)) UInt p, UInt dims[1]) {
dims[0] = m;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<2> {
static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p,
UInt dims[2]) {
dims[0] = m;
dims[1] = n;
}
};
/* -------------------------------------------------------------------------- */
template <> struct DimHelper<3> {
static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) {
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
/**
* @class TensorProxy aka_types.hh
* @desc 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 {
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;
}
#ifndef SWIG
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);
}
#endif
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; }
#ifndef SWIG
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;
}
#endif
// 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;
UInt n[ndim];
};
/* -------------------------------------------------------------------------- */
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;
// }
/* ------------------------------------------------------------------------ */
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);
this->_size = src.size();
}
TensorStorage() : values(nullptr) {
for (UInt d = 0; d < ndim; ++d)
this->n[d] = 0;
_size = 0;
}
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];
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];
static_assert(std::is_trivially_copyable<T>{},
"Cannot copy a tensor on non trivial types");
memcpy((void *)this->values, (void *)src.storage(),
this->_size * sizeof(T));
this->wrapped = false;
}
virtual ~TensorStorage() {
if (!this->wrapped)
delete[] this->values;
}
/* ------------------------------------------------------------------------ */
inline TensorStorage & operator=(const TensorStorage & src) {
return this->operator=(dynamic_cast<RetType &>(src));
}
/* ------------------------------------------------------------------------ */
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");
memcpy((void *)this->values, (void *)src.storage(),
this->_size * sizeof(T));
} 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));
}
/// Y = \alpha X + Y
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 clear() { memset(values, 0, _size * sizeof(T)); };
inline void set(const T & t) { std::fill_n(values, _size, t); };
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");
memcpy(values, other.storage(), _size * sizeof(T));
}
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);
}
};
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:
UInt n[ndim];
UInt _size;
T * values;
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:
~Vector() override = default;
/* ------------------------------------------------------------------------ */
inline Vector & operator=(const Vector & src) {
parent::operator=(src);
return *this;
}
/* ------------------------------------------------------------------------ */
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) {
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, Real alpha = 1.0);
/* ------------------------------------------------------------------------ */
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)
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; }
/* ------------------------------------------------------------------------ */
/// 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 << "]";
}
// friend class ::akantu::Array<T>;
};
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;
}
/* ------------------------------------------------------------------------ */
/* 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, 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;
}
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 const 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]);
}
- inline void block(Matrix & block, UInt pos_i, UInt pos_j) {
+ 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) {
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) 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());
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->clear();
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);
else if (this->cols() == 2)
return Math::det2(this->values);
else if (this->cols() == 3)
return Math::det3(this->values);
else
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());
else if (this->cols() == 2)
return Math::matrixDoubleDot22(this->values, other.storage());
else if (this->cols() == 3)
return Math::matrixDoubleDot33(this->values, other.storage());
else
AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions"
<< " than 1, 2 or 3.");
return T();
}
/* ---------------------------------------------------------------------- */
/// 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,
Real 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 const 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 const 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
#endif /* __AKANTU_AKA_TYPES_HH__ */
diff --git a/src/fe_engine/element_class_structural.hh b/src/fe_engine/element_class_structural.hh
index 1e3a79b47..3dbefe93d 100644
--- a/src/fe_engine/element_class_structural.hh
+++ b/src/fe_engine/element_class_structural.hh
@@ -1,277 +1,250 @@
/**
* @file element_class_structural.hh
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Damien Spielmann <damien.spielmann@epfl.ch>
*
* @date creation: Thu Feb 21 2013
* @date last modification: Thu Oct 22 2015
*
* @brief Specialization of the element classes for structural elements
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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 "element_class.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__
#define __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__
namespace akantu {
/// Macro to generate the InterpolationProperty structures for different
/// interpolation types
#define AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY( \
- itp_type, itp_geom_type, ndof, nb_stress) \
+ itp_type, itp_geom_type, ndof, nb_stress, nb_dnds_cols) \
template <> struct InterpolationProperty<itp_type> { \
static const InterpolationKind kind{_itk_structural}; \
static const UInt nb_nodes_per_element{ \
InterpolationProperty<itp_geom_type>::nb_nodes_per_element}; \
static const InterpolationType itp_geometry_type{itp_geom_type}; \
static const UInt natural_space_dimension{ \
InterpolationProperty<itp_geom_type>::natural_space_dimension}; \
static const UInt nb_degree_of_freedom{ndof}; \
static const UInt nb_stress_components{nb_stress}; \
+ static const UInt dnds_columns{nb_dnds_cols}; \
}
/* -------------------------------------------------------------------------- */
template <InterpolationType interpolation_type>
class InterpolationElement<interpolation_type, _itk_structural> {
public:
using interpolation_property = InterpolationProperty<interpolation_type>;
/// compute the shape values for a given set of points in natural coordinates
static inline void computeShapes(const Matrix<Real> & natural_coord,
const Matrix<Real> & real_coord,
Tensor3<Real> & N) {
for (UInt i = 0; i < natural_coord.cols(); ++i) {
Matrix<Real> n_t = N(i);
computeShapes(natural_coord(i), real_coord, n_t);
}
}
/// compute the shape values for a given point in natural coordinates
static inline void computeShapes(const Vector<Real> & natural_coord,
const Matrix<Real> & real_coord,
Matrix<Real> & N);
/// compute shape derivatives (input is dxds) for a set of points
static inline void computeShapeDerivatives(const Tensor3<Real> & Js,
const Tensor3<Real> & DNDSs,
const Matrix<Real> & R,
Tensor3<Real> & Bs) {
for (UInt i = 0; i < Js.size(2); ++i) {
Matrix<Real> J = Js(i);
Matrix<Real> DNDS = DNDSs(i);
- Matrix<Real> DNDS_R(DNDS.rows(), DNDS.cols());
- DNDS_R.mul<false, false>(DNDS, R);
- Matrix<Real> B = Bs(i);
+ Matrix<Real> DNDX(DNDS.rows(), DNDS.cols());
auto inv_J = J.inverse();
- Matrix<Real> inv_J_full(DNDS.rows(), DNDS.rows());
-
- // Gotta repeat J^-1 for each stress component
- for (UInt k = 0, pos = 0; k < getNbStressComponents();
- k++, pos += inv_J.rows()) {
- inv_J_full.block(inv_J, pos, pos);
- }
-
- B.mul<false, false>(inv_J_full, DNDS_R);
+ DNDX.mul<false, false>(inv_J, DNDS);
+ Matrix<Real> B_R = Bs(i);
+ Matrix<Real> B(B_R.rows(), B_R.cols());
+ arrangeInVoigt(DNDX, B);
+ B_R.mul<false, false>(B, R);
}
}
/**
* compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of
* shape functions along with variation of natural coordinates on a given set
* of points in natural coordinates
*/
static inline void computeDNDS(const Matrix<Real> & natural_coord,
const Matrix<Real> & real_coord,
Tensor3<Real> & dnds) {
for (UInt i = 0; i < natural_coord.cols(); ++i) {
Matrix<Real> dnds_t = dnds(i);
computeDNDS(natural_coord(i), real_coord, dnds_t);
}
}
/**
* compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of
* shape functions along with
* variation of natural coordinates on a given point in natural
* coordinates
*/
static inline void computeDNDS(const Vector<Real> & natural_coord,
const Matrix<Real> & real_coord,
Matrix<Real> & dnds);
+ /**
+ * arrange B in Voigt notation from DNDS
+ */
+ static inline void arrangeInVoigt(const Matrix<Real> & dnds,
+ Matrix<Real> & B) {
+ // Default implementation assumes dnds is already in Voigt notation
+ B.deepCopy(dnds);
+ }
+
public:
static AKANTU_GET_MACRO_NOT_CONST(
NbNodesPerInterpolationElement,
interpolation_property::nb_nodes_per_element, UInt);
static AKANTU_GET_MACRO_NOT_CONST(
ShapeSize, (interpolation_property::nb_nodes_per_element *
interpolation_property::nb_degree_of_freedom *
interpolation_property::nb_degree_of_freedom),
UInt);
static AKANTU_GET_MACRO_NOT_CONST(
ShapeDerivativesSize, (interpolation_property::nb_nodes_per_element *
interpolation_property::nb_degree_of_freedom *
interpolation_property::nb_stress_components),
UInt);
static AKANTU_GET_MACRO_NOT_CONST(
NaturalSpaceDimension, interpolation_property::natural_space_dimension,
UInt);
static AKANTU_GET_MACRO_NOT_CONST(
NbDegreeOfFreedom, interpolation_property::nb_degree_of_freedom, UInt);
static AKANTU_GET_MACRO_NOT_CONST(
NbStressComponents, interpolation_property::nb_stress_components, UInt);
};
/// Macro to generate the element class structures for different structural
/// element types
/* -------------------------------------------------------------------------- */
#define AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY( \
elem_type, geom_type, interp_type, parent_el_type, elem_kind, sp, \
gauss_int_type, min_int_order) \
template <> struct ElementClassProperty<elem_type> { \
static const GeometricalType geometrical_type{geom_type}; \
static const InterpolationType interpolation_type{interp_type}; \
static const ElementType parent_element_type{parent_el_type}; \
static const ElementKind element_kind{elem_kind}; \
static const UInt spatial_dimension{sp}; \
static const GaussIntegrationType gauss_integration_type{gauss_int_type}; \
static const UInt polynomial_degree{min_int_order}; \
}
/* -------------------------------------------------------------------------- */
/* ElementClass for structural elements */
/* -------------------------------------------------------------------------- */
template <ElementType element_type>
class ElementClass<element_type, _ek_structural>
: public GeometricalElement<
ElementClassProperty<element_type>::geometrical_type>,
public InterpolationElement<
ElementClassProperty<element_type>::interpolation_type> {
protected:
using geometrical_element =
GeometricalElement<ElementClassProperty<element_type>::geometrical_type>;
using interpolation_element = InterpolationElement<
ElementClassProperty<element_type>::interpolation_type>;
using parent_element =
ElementClass<ElementClassProperty<element_type>::parent_element_type>;
public:
static inline void
computeRotationMatrix(Matrix<Real> & /*R*/, const Matrix<Real> & /*X*/,
const Vector<Real> & /*extra_normal*/) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute jacobian (or integration variable change factor) for a given point
- static inline void computeJMat(const Matrix<Real> & DNDS,
+ static inline void computeJMat(const Vector<Real> & natural_coords,
const Matrix<Real> & Xs, Matrix<Real> & J) {
- auto nb_nodes = Xs.cols();
- auto dim = Xs.rows();
- auto nb_dof =
- interpolation_element::interpolation_property::nb_degree_of_freedom;
- Matrix<Real> B(dim, nb_nodes);
- for (UInt s = 0; s < dim; ++s) {
- for (UInt n = 0; n < nb_nodes; ++n) {
- B(s, n) = DNDS(s, n * nb_dof + s);
- }
- }
-
- J.mul<false, true>(B, Xs);
+ Matrix<Real> dnds(Xs.rows(), Xs.cols());
+ parent_element::computeDNDS(natural_coords, dnds);
+ J.mul<false, true>(dnds, Xs);
}
- static inline void computeJMat(const Tensor3<Real> & DNDSs,
+ static inline void computeJMat(const Matrix<Real> & natural_coords,
const Matrix<Real> & Xs, Tensor3<Real> & Js) {
- using itp = typename interpolation_element::interpolation_property;
- auto nb_nodes = Xs.cols();
- auto dim = Xs.rows();
- auto nb_dof = itp::nb_degree_of_freedom;
- Matrix<Real> B(dim, nb_nodes);
-
- for (UInt i = 0; i < Xs.cols(); ++i) {
- Matrix<Real> DNDS = DNDSs(i);
+ for (UInt i = 0 ; i < natural_coords.cols(); ++i) {
+ // because non-const l-value reference does not bind to r-value
Matrix<Real> J = Js(i);
-
- B.clear();
- for (UInt s = 0; s < dim; ++s) {
- for (UInt n = 0; n < nb_nodes; ++n) {
- B(s, n) = DNDS(s, n * nb_dof + s);
- }
- }
-
- parent_element::computeJMat(B, Xs, J);
- // computeJMat(DNDS, Xs, J);
+ computeJMat(Vector<Real>(natural_coords(i)), Xs, J);
}
}
static inline void computeJacobian(const Matrix<Real> & natural_coords,
const Matrix<Real> & node_coords,
Vector<Real> & jacobians) {
using itp = typename interpolation_element::interpolation_property;
- UInt nb_points = natural_coords.cols();
- Matrix<Real> dnds(itp::natural_space_dimension, itp::nb_nodes_per_element);
- Matrix<Real> J(natural_coords.rows(), itp::natural_space_dimension);
-
- // Extract relevant first lines
- auto x = node_coords.block(0, 0, itp::natural_space_dimension,
- itp::nb_nodes_per_element);
-
- for (UInt p = 0; p < nb_points; ++p) {
- Vector<Real> ncoord_p(natural_coords(p));
- parent_element::computeDNDS(ncoord_p, dnds);
- parent_element::computeJMat(dnds, x, J);
- jacobians(p) = J.det();
+ Tensor3<Real> Js(itp::natural_space_dimension, itp::natural_space_dimension,
+ natural_coords.cols());
+ computeJMat(natural_coords, node_coords, Js);
+ for (UInt i = 0; i < natural_coords.cols(); ++i) {
+ Matrix<Real> J = Js(i);
+ jacobians(i) = J.det();
}
}
static inline void computeRotation(const Matrix<Real> & node_coords,
Matrix<Real> & rotation);
public:
static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_structural, ElementKind);
static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, _not_defined, ElementType);
static AKANTU_GET_MACRO_NOT_CONST(FacetType, _not_defined, ElementType);
static constexpr auto getFacetType(__attribute__((unused)) UInt t = 0) {
return _not_defined;
}
static constexpr AKANTU_GET_MACRO_NOT_CONST(
SpatialDimension, ElementClassProperty<element_type>::spatial_dimension,
UInt);
static constexpr auto getFacetTypes() {
return ElementClass<_not_defined>::getFacetTypes();
}
};
} // namespace akantu
/* -------------------------------------------------------------------------- */
#include "element_classes/element_class_hermite_inline_impl.cc"
/* keep order */
#include "element_classes/element_class_bernoulli_beam_inline_impl.cc"
#include "element_classes/element_class_kirchhoff_shell_inline_impl.cc"
/* -------------------------------------------------------------------------- */
#endif /* __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__ */
diff --git a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc
index f7ca15021..407cb0af0 100644
--- a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc
+++ b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc
@@ -1,202 +1,222 @@
/**
* @file element_class_bernoulli_beam_inline_impl.cc
*
* @author Fabian Barras <fabian.barras@epfl.ch>
*
* @date creation: Fri Jul 15 2011
* @date last modification: Sun Oct 19 2014
*
* @brief Specialization of the element_class class for the type
_bernoulli_beam_2
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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/>.
*
* @section DESCRIPTION
*
* @verbatim
--x-----q1----|----q2-----x---> x
-1 0 1
@endverbatim
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_static_if.hh"
#include "element_class_structural.hh"
//#include "aka_element_classes_info.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__
#define __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_2,
_itp_lagrange_segment_2, 3,
- 2);
+ 2, 6);
AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_3,
_itp_lagrange_segment_2, 6,
- 4);
+ 4, 6);
AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_2,
_gt_segment_2,
_itp_bernoulli_beam_2,
_segment_2, _ek_structural, 2,
_git_segment, 3);
AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_3,
_gt_segment_2,
_itp_bernoulli_beam_3,
_segment_2, _ek_structural, 3,
_git_segment, 3);
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeShapes(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
Matrix<Real> & N) {
Vector<Real> L(2);
InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes(
natural_coords, L);
Matrix<Real> H(2, 4);
InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes(
natural_coords, real_coord, H);
// clang-format off
// u1 v1 t1 u2 v2 t2
N = {{L(0), 0 , 0 , L(1), 0 , 0 }, // u
{0 , H(0, 0), H(0, 1), 0 , H(0, 2), H(0, 3)}, // v
{0 , H(1, 0), H(1, 1), 0 , H(1, 2), H(1, 3)}}; // theta
// clang-format on
}
template <>
inline void
InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeShapes(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
Matrix<Real> & N) {
Vector<Real> L(2);
InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes(
natural_coords, L);
Matrix<Real> H(2, 4);
InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes(
natural_coords, real_coord, H);
// clang-format off
// u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2
N = {{L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 , 0 , 0 , 0 }, // u
{0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3), 0 }, // v
{0 , 0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3)}, // w
{0 , 0 , 0 , L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 }, // thetax
{0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3), 0 }, // thetay
{0 , 0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3)}}; // thetaz
// clang-format on
}
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeDNDS(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
- Matrix<Real> & B) {
+ Matrix<Real> & dnds) {
Matrix<Real> L(1, 2);
InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeDNDS(
natural_coords, L);
Matrix<Real> H(1, 4);
InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS(
natural_coords, real_coord, H);
+ // Storing the derivatives in dnds
+ dnds.block(L, 0, 0);
+ dnds.block(H, 0, 2);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+inline void
+InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::arrangeInVoigt(
+ const Matrix<Real> & dnds, Matrix<Real> & B) {
+ auto L = dnds.block(0, 0, 1, 2); // Lagrange shape derivatives
+ auto H = dnds.block(0, 2, 1, 4); // Hermite shape derivatives
// clang-format off
- // u1 v1 t1 u2 v2 t2
- B = {{L(0, 0), 0 , 0 , L(0, 1), 0 , 0 }, // epsilon
- {0 , H(0, 0), H(0, 1), 0 , H(0, 2), H(0, 3)}}; // chi (curv.)
+ // u1 v1 t1 u2 v2 t2
+ B = {{L(0, 0), 0, 0, L(0, 1), 0, 0 },
+ {0, H(0, 0), H(0, 1), 0, H(0, 2), H(0, 3)}};
// clang-format on
}
+/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeDNDS(
- const Vector<Real> & natural_coords, const Matrix<Real> & real_coord, Matrix<Real> & B) {
- Matrix<Real> L(1, 2);
- InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeDNDS(
- natural_coords, L);
- Matrix<Real> H(1, 4);
- InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS(
- natural_coords, real_coord, H);
+ const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
+ Matrix<Real> & dnds) {
+ InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeDNDS(
+ natural_coords, real_coord, dnds);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+inline void
+InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::arrangeInVoigt(
+ const Matrix<Real> & dnds, Matrix<Real> & B) {
+ auto L = dnds.block(0, 0, 1, 2); // Lagrange shape derivatives
+ auto H = dnds.block(0, 2, 1, 4); // Hermite shape derivatives
+
// clang-format off
// u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2
B = {{L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 , 0 , 0 , 0 }, // eps
{0 , H(0, 0), 0 , 0 , 0 , H(0, 1), 0 , H(0, 2), 0 , 0 , 0 , H(0, 3)}, // chi strong axis
{0 , 0 ,-H(0, 0), 0 , H(0, 1), 0 , 0 , 0 ,-H(0, 2), 0 , H(0, 3), 0 }, // chi weak axis
{0 , 0 , 0 , L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 }}; // chi torsion
// clang-format on
}
/* -------------------------------------------------------------------------- */
template <>
inline void ElementClass<_bernoulli_beam_2>::computeRotationMatrix(
Matrix<Real> & R, const Matrix<Real> & X, const Vector<Real> &) {
Vector<Real> x2 = X(1); // X2
Vector<Real> x1 = X(0); // X1
auto cs = (x2 - x1);
cs.normalize();
auto c = cs(0);
auto s = cs(1);
// clang-format off
/// Definition of the rotation matrix
R = {{ c, s, 0.},
{-s, c, 0.},
{ 0., 0., 1.}};
// clang-format on
}
/* -------------------------------------------------------------------------- */
template <>
inline void ElementClass<_bernoulli_beam_3>::computeRotationMatrix(
Matrix<Real> & R, const Matrix<Real> & X, const Vector<Real> & n) {
Vector<Real> x2 = X(1); // X2
Vector<Real> x1 = X(0); // X1
auto dim = X.rows();
auto x = (x2 - x1);
x.normalize();
auto x_n = x.crossProduct(n);
Matrix<Real> Pe = {{1., 0., 0.}, {0., -1., 0.}, {0., 0., 1.}};
Matrix<Real> Pg(dim, dim);
Pg(0) = x;
Pg(1) = x_n;
Pg(2) = n;
Pe *= Pg.inverse();
R.clear();
/// Definition of the rotation matrix
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
R(i + dim, j + dim) = R(i, j) = Pe(i, j);
}
} // namespace akantu
#endif /* __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__ */
diff --git a/src/fe_engine/element_classes/element_class_hermite_inline_impl.cc b/src/fe_engine/element_classes/element_class_hermite_inline_impl.cc
index ba268905d..8d79e01c5 100644
--- a/src/fe_engine/element_classes/element_class_hermite_inline_impl.cc
+++ b/src/fe_engine/element_classes/element_class_hermite_inline_impl.cc
@@ -1,180 +1,180 @@
/**
* @file element_class_hermite_inline_impl.cc
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Fri Jul 15 2011
* @date last modification: Sun Oct 19 2014
*
* @brief Specialization of the element_class class for the type
_hermite
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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/>.
*
* @section DESCRIPTION
*
* @verbatim
--x-----q1----|----q2-----x---> x
-1 0 1
@endverbatim
*
* @subsection shapes Shape functions
* @f[
* \begin{array}{ll}
* M_1(\xi) &= 1/4(\xi^{3}/-3\xi+2)\\
* M_2(\xi) &= -1/4(\xi^{3}-3\xi-2)
* \end{array}
*
* \begin{array}{ll}
* L_1(\xi) &= 1/4(\xi^{3}-\xi^{2}-\xi+1)\\
* L_2(\xi) &= 1/4(\xi^{3}+\xi^{2}-\xi-1)
* \end{array}
*
* \begin{array}{ll}
* M'_1(\xi) &= 3/4(\xi^{2}-1)\\
* M'_2(\xi) &= -3/4(\xi^{2}-1)
* \end{array}
*
* \begin{array}{ll}
* L'_1(\xi) &= 1/4(3\xi^{2}-2\xi-1)\\
* L'_2(\xi) &= 1/4(3\xi^{2}+2\xi-1)
* \end{array}
*@f]
*
* @subsection dnds Shape derivatives
*
*@f[
* \begin{array}{ll}
* N'_1(\xi) &= -1/2\\
* N'_2(\xi) &= 1/2
* \end{array}]
*
* \begin{array}{ll}
* -M''_1(\xi) &= -3\xi/2\\
* -M''_2(\xi) &= 3\xi/2\\
* \end{array}
*
* \begin{array}{ll}
* -L''_1(\xi) &= -1/2a(3\xi/a-1)\\
* -L''_2(\xi) &= -1/2a(3\xi/a+1)
* \end{array}
*@f]
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_static_if.hh"
#include "element_class_structural.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_CC__
#define __AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_hermite_2,
_itp_lagrange_segment_2, 2,
- 1);
+ 1, 4);
/* -------------------------------------------------------------------------- */
namespace {
namespace details {
inline Real computeLength(const Matrix<Real> & real_coord) {
Vector<Real> x1 = real_coord(0);
Vector<Real> x2 = real_coord(1);
return x1.distance(x2);
}
inline void computeShapes(const Vector<Real> & natural_coords, Real a, Matrix<Real> & N) {
/// natural coordinate
Real xi = natural_coords(0);
// Cubic Hermite splines interpolating displacement
auto M1 = 1. / 4. * Math::pow<2>(xi - 1) * (xi + 2);
auto M2 = 1. / 4. * Math::pow<2>(xi + 1) * (2 - xi);
auto L1 = a / 4. * Math::pow<2>(xi - 1) * (xi + 1);
auto L2 = a / 4. * Math::pow<2>(xi + 1) * (xi - 1);
#if 1 // Version where we also interpolate the rotations
// Derivatives (with respect to x) of previous functions interpolating
// rotations
auto M1_ = 3. / (4. * a) * (xi * xi - 1);
auto M2_ = 3. / (4. * a) * (1 - xi * xi);
auto L1_ = 1 / 4. *
(3 * xi * xi - 2 * xi - 1);
auto L2_ = 1 / 4. *
(3 * xi * xi + 2 * xi - 1);
// clang-format off
// v1 t1 v2 t2
N = {{M1 , L1 , M2 , L2}, // displacement interpolation
{M1_, L1_, M2_, L2_}}; // rotation interpolation
// clang-format on
#else // Version where we only interpolate displacements
// clang-format off
// v1 t1 v2 t2
N = {{M1, L1, M2, L2}};
// clang-format on
#endif
}
/* ---------------------------------------------------------------------- */
inline void computeDNDS(const Vector<Real> & natural_coords, Real a,
Matrix<Real> & B) {
// natural coordinate
Real xi = natural_coords(0);
// Derivatives with respect to xi for rotations
auto M1__ = 3. / 2. * xi;
auto M2__ = 3. / 2. * (-xi);
auto L1__ = a / 2. * (3 * xi - 1);
auto L2__ = a / 2. * (3 * xi + 1);
// v1 t1 v2 t2
B = {{M1__, L1__, M2__, L2__}}; // computing curvature : {chi} = [B]{d}
B /= a; // to account for first order deriv w/r to x
}
} // namespace details
} // namespace
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
Matrix<Real> & N) {
auto L = details::computeLength(real_coord);
details::computeShapes(natural_coords, L / 2, N);
}
/* -------------------------------------------------------------------------- */
template <>
inline void InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coord,
Matrix<Real> & B) {
auto L = details::computeLength(real_coord);
details::computeDNDS(natural_coords, L / 2, B);
}
} // namespace akantu
#endif /* __AKANTU_ELEMENT_CLASS_HERMITE_INLINE_IMPL_CC__ */
diff --git a/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc b/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc
index 65b452d57..6818c8ddc 100644
--- a/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc
+++ b/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc
@@ -1,189 +1,211 @@
/**
* @file element_class_kirchhoff_shell_inline_impl.cc
*
* @author Damien Spielmann <damien.spielmann@epfl.ch>
*
* @date creation: Fri Jul 04 2014
* @date last modification: Sun Oct 19 2014
*
* @brief Element class Kirchhoff Shell
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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 "element_class_structural.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__
#define __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(
- _itp_discrete_kirchhoff_triangle_18, _itp_lagrange_triangle_3, 6, 6);
+ _itp_discrete_kirchhoff_triangle_18, _itp_lagrange_triangle_3, 6, 6, 21);
AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(
_discrete_kirchhoff_triangle_18, _gt_triangle_3,
_itp_discrete_kirchhoff_triangle_18, _triangle_3, _ek_structural, 3,
_git_triangle, 2);
/* -------------------------------------------------------------------------- */
namespace detail {
inline void computeBasisChangeMatrix(Matrix<Real> & P,
const Matrix<Real> & X) {
Vector<Real> X1 = X(0);
Vector<Real> X2 = X(1);
Vector<Real> X3 = X(2);
Vector<Real> a1 = X2 - X1;
Vector<Real> a2 = X3 - X1;
a1.normalize();
Vector<Real> e3 = a1.crossProduct(a2);
e3.normalize();
Vector<Real> e2 = e3.crossProduct(a1);
P(0) = a1;
P(1) = e2;
P(2) = e3;
}
}
/* -------------------------------------------------------------------------- */
template <>
inline void
ElementClass<_discrete_kirchhoff_triangle_18>::computeRotationMatrix(
Matrix<Real> & R, const Matrix<Real> & X, const Vector<Real> &) {
auto dim = X.rows();
Matrix<Real> P(dim, dim);
detail::computeBasisChangeMatrix(P, X);
R.clear();
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
R(i + dim, j + dim) = R(i, j) = P(j, i);
}
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_discrete_kirchhoff_triangle_18>::computeShapes(
const Vector<Real> & /*natural_coords*/,
const Matrix<Real> & /*real_coord*/, Matrix<Real> & /*N*/) {}
/* -------------------------------------------------------------------------- */
template <>
inline void
InterpolationElement<_itp_discrete_kirchhoff_triangle_18>::computeDNDS(
const Vector<Real> & natural_coords, const Matrix<Real> & real_coordinates,
Matrix<Real> & B) {
auto dim = real_coordinates.cols();
Matrix<Real> P(dim, dim);
detail::computeBasisChangeMatrix(P, real_coordinates);
auto X = P * real_coordinates;
Vector<Real> X1 = X(0);
Vector<Real> X2 = X(1);
Vector<Real> X3 = X(2);
std::array<Vector<Real>, 3> A = {X2 - X1, X3 - X2, X1 - X3};
std::array<Real, 3> L, C, S;
// Setting all last coordinates to 0
std::for_each(A.begin(), A.end(), [](auto & a) { a(2) = 0; });
// Computing lengths
std::transform(A.begin(), A.end(), L.begin(),
[](Vector<Real> & a) { return a.norm<L_2>(); });
// Computing cosines
std::transform(A.begin(), A.end(), L.begin(), C.begin(),
[](Vector<Real> & a, Real & l) { return a(0) / l; });
// Computing sines
std::transform(A.begin(), A.end(), L.begin(), S.begin(),
[](Vector<Real> & a, Real & l) { return a(1) / l; });
// Natural coordinates
Real xi = natural_coords(0);
Real eta = natural_coords(1);
// Derivative of quadratic interpolation functions
Matrix<Real> dP = {{4 * (1 - 2 * xi - eta), 4 * eta, -4 * eta},
{-4 * xi, 4 * xi, 4 * (1 - xi - 2 * eta)}};
Matrix<Real> dNx1 = {
{3. / 2 * (dP(0, 0) * C[0] / L[0] - dP(0, 2) * C[2] / L[2]),
3. / 2 * (dP(0, 1) * C[1] / L[1] - dP(0, 0) * C[0] / L[0]),
3. / 2 * (dP(0, 2) * C[2] / L[2] - dP(0, 1) * C[1] / L[1])},
{3. / 2 * (dP(1, 0) * C[0] / L[0] - dP(1, 2) * C[2] / L[2]),
3. / 2 * (dP(1, 1) * C[1] / L[1] - dP(1, 0) * C[0] / L[0]),
3. / 2 * (dP(1, 2) * C[2] / L[2] - dP(1, 1) * C[1] / L[1])}};
Matrix<Real> dNx2 = {
// clang-format off
{-1 - 3. / 4 * (dP(0, 0) * C[0] * C[0] + dP(0, 2) * C[2] * C[2]),
1 - 3. / 4 * (dP(0, 1) * C[1] * C[1] + dP(0, 0) * C[0] * C[0]),
-3. / 4 * (dP(0, 2) * C[2] * C[2] + dP(0, 1) * C[1] * C[1])},
{-1 - 3. / 4 * (dP(1, 0) * C[0] * C[0] + dP(1, 2) * C[2] * C[2]),
-3. / 4 * (dP(1, 1) * C[1] * C[1] + dP(1, 0) * C[0] * C[0]),
1 - 3. / 4 * (dP(1, 2) * C[2] * C[2] + dP(1, 1) * C[1] * C[1])}};
// clang-format on
Matrix<Real> dNx3 = {
{-3. / 4 * (dP(0, 0) * C[0] * S[0] + dP(0, 2) * C[2] * S[2]),
-3. / 4 * (dP(0, 1) * C[1] * S[1] + dP(0, 0) * C[0] * S[0]),
-3. / 4 * (dP(0, 2) * C[2] * S[2] + dP(0, 1) * C[1] * S[1])},
{-3. / 4 * (dP(1, 0) * C[0] * S[0] + dP(1, 2) * C[2] * S[2]),
-3. / 4 * (dP(1, 1) * C[1] * S[1] + dP(1, 0) * C[0] * S[0]),
-3. / 4 * (dP(1, 2) * C[2] * S[2] + dP(1, 1) * C[1] * S[1])}};
Matrix<Real> dNy1 = {
{3. / 2 * (dP(0, 0) * S[0] / L[0] - dP(0, 2) * S[2] / L[2]),
3. / 2 * (dP(0, 1) * S[1] / L[1] - dP(0, 0) * S[0] / L[0]),
3. / 2 * (dP(0, 2) * S[2] / L[2] - dP(0, 1) * S[1] / L[1])},
{3. / 2 * (dP(1, 0) * S[0] / L[0] - dP(1, 2) * S[2] / L[2]),
3. / 2 * (dP(1, 1) * S[1] / L[1] - dP(1, 0) * S[0] / L[0]),
3. / 2 * (dP(1, 2) * S[2] / L[2] - dP(1, 1) * S[1] / L[1])}};
Matrix<Real> dNy2 = dNx3;
Matrix<Real> dNy3 = {
// clang-format off
{-1 - 3. / 4 * (dP(0, 0) * S[0] * S[0] + dP(0, 2) * S[2] * S[2]),
1 - 3. / 4 * (dP(0, 1) * S[1] * S[1] + dP(0, 0) * S[0] * S[0]),
-3. / 4 * (dP(0, 2) * S[2] * S[2] + dP(0, 1) * S[1] * S[1])},
{-1 - 3. / 4 * (dP(1, 0) * S[0] * S[0] + dP(1, 2) * S[2] * S[2]),
-3. / 4 * (dP(1, 1) * S[1] * S[1] + dP(1, 0) * S[0] * S[0]),
1 - 3. / 4 * (dP(1, 2) * S[2] * S[2] + dP(1, 1) * S[1] * S[1])}};
// clang-format on
// Derivative of linear (membrane mode) functions
- Matrix<Real> dNm = {{1, -1, 0}, {-1, 0, 1}};
+ Matrix<Real> dNm(2, 3);
+ InterpolationElement<_itp_lagrange_triangle_3, _itk_lagrangian>::computeDNDS(
+ natural_coords, dNm);
+
+ UInt i = 0;
+ for (const Matrix<Real> & mat : {dNm, dNx1, dNx2, dNx3, dNy1, dNy2, dNy3}) {
+ B.block(mat, 0, i);
+ i += mat.cols();
+ }
+}
+/* -------------------------------------------------------------------------- */
+template <>
+inline void
+InterpolationElement<_itp_discrete_kirchhoff_triangle_18,
+ _itk_structural>::arrangeInVoigt(const Matrix<Real> & dnds,
+ Matrix<Real> & B) {
+ Matrix<Real> dNm(2, 3), dNx1(2, 3), dNx2(2, 3), dNx3(2, 3), dNy1(2, 3),
+ dNy2(2, 3), dNy3(2, 3);
+ UInt i = 0;
+ for (Matrix<Real> * mat : {&dNm, &dNx1, &dNx2, &dNx3, &dNy1, &dNy2, &dNy3}) {
+ *mat = dnds.block(0, i, 2, 3);
+ i += mat->cols();
+ }
// clang-format off
B = {
// Membrane shape functions; TODO use the triangle_3 shapes
{dNm(0, 0), 0, 0, 0, 0, 0, dNm(0, 1), 0, 0, 0, 0, 0, dNm(0, 2), 0, 0, 0, 0, 0, },
{0, dNm(1, 0), 0, 0, 0, 0, 0, dNm(1, 1), 0, 0, 0, 0, 0, dNm(1, 2), 0, 0, 0, 0, },
{dNm(1, 0), dNm(0, 0), 0, 0, 0, 0, dNm(1, 1), dNm(0, 1), 0, 0, 0, 0, dNm(1, 2), dNm(0, 2), 0, 0, 0, 0, },
// Bending shape functions
{0, 0, dNx1(0, 0), -dNx3(0, 0), dNx2(0, 0), 0, 0, 0, dNx1(0, 1), -dNx3(0, 1), dNx2(0, 1), 0, 0, 0, dNx1(0, 2), -dNx3(0, 2), dNx2(0, 2), 0},
{0, 0, dNy1(1, 0), -dNy3(1, 0), dNy2(1, 0), 0, 0, 0, dNy1(1, 1), -dNy3(1, 1), dNy2(1, 1), 0, 0, 0, dNy1(1, 2), -dNy3(1, 2), dNy2(1, 2), 0},
{0, 0, dNx1(1, 0) + dNy1(0, 0), -dNx3(1, 0) - dNy3(0, 0), dNx2(1, 0) + dNy2(1, 0), 0, 0, 0, dNx1(1, 1) + dNy1(0, 1), -dNx3(1, 1) - dNy3(0, 0), dNx2(1, 1) + dNy2(1, 1), 0, 0, 0, dNx1(1, 2) + dNy1(0, 2), -dNx3(1, 2) - dNy3(0, 2), dNx2(1, 2) + dNy2(1, 2), 0}};
// clang-format on
}
} // namespace akantu
#endif /* __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__ */
diff --git a/src/fe_engine/shape_structural_inline_impl.cc b/src/fe_engine/shape_structural_inline_impl.cc
index 0eefa86d8..5792b0d34 100644
--- a/src/fe_engine/shape_structural_inline_impl.cc
+++ b/src/fe_engine/shape_structural_inline_impl.cc
@@ -1,383 +1,389 @@
/**
* @file shape_structural_inline_impl.cc
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Dec 13 2010
* @date last modification: Thu Oct 15 2015
*
* @brief ShapeStructural inline implementation
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "mesh_iterators.hh"
#include "shape_structural.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__
#define __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__
namespace akantu {
namespace {
/// Extract nodal coordinates per elements
template <ElementType type>
std::unique_ptr<Array<Real>>
getNodesPerElement(const Mesh & mesh, const Array<Real> & nodes,
const GhostType & ghost_type) {
const auto dim = ElementClass<type>::getSpatialDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
auto nodes_per_element = std::make_unique<Array<Real>>(0, dim * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, *nodes_per_element, type, ghost_type);
return nodes_per_element;
}
}
template <ElementKind kind>
inline void ShapeStructural<kind>::initShapeFunctions(
const Array<Real> & /* unused */, const Matrix<Real> & /* unused */,
const ElementType & /* unused */, const GhostType & /* unused */) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
#define INIT_SHAPE_FUNCTIONS(type) \
setIntegrationPointsByType<type>(integration_points, ghost_type); \
precomputeRotationMatrices<type>(nodes, ghost_type); \
precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type); \
precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);
template <>
inline void ShapeStructural<_ek_structural>::initShapeFunctions(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
const ElementType & type, const GhostType & ghost_type) {
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
}
#undef INIT_SHAPE_FUNCTIONS
/* -------------------------------------------------------------------------- */
template <>
template <ElementType type>
void ShapeStructural<_ek_structural>::computeShapesOnIntegrationPoints(
const Array<Real> & nodes, const Matrix<Real> & integration_points,
Array<Real> & shapes, const GhostType & ghost_type,
const Array<UInt> & filter_elements) const {
UInt nb_points = integration_points.cols();
UInt nb_element = mesh.getConnectivity(type, ghost_type).size();
shapes.resize(nb_element * nb_points);
UInt ndof = ElementClass<type>::getNbDegreeOfFreedom();
#if !defined(AKANTU_NDEBUG)
UInt size_of_shapes = ElementClass<type>::getShapeSize();
AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes,
"The shapes array does not have the correct "
<< "number of component");
#endif
auto shapes_it = shapes.begin_reinterpret(
ElementClass<type>::getNbNodesPerInterpolationElement(), ndof, nb_points,
nb_element);
auto shapes_begin = shapes_it;
if (filter_elements != empty_filter) {
nb_element = filter_elements.size();
}
auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
auto nodes_it = nodes_per_element->begin(mesh.getSpatialDimension(),
Mesh::getNbNodesPerElement(type));
auto nodes_begin = nodes_it;
for (UInt elem = 0; elem < nb_element; ++elem) {
if (filter_elements != empty_filter) {
shapes_it = shapes_begin + filter_elements(elem);
nodes_it = nodes_begin + filter_elements(elem);
}
Tensor3<Real> & N = *shapes_it;
auto & real_coord = *nodes_it;
ElementClass<type>::computeShapes(integration_points, real_coord, N);
if (filter_elements == empty_filter)
++shapes_it;
}
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeRotationMatrices(
const Array<Real> & nodes, const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const auto spatial_dimension = mesh.getSpatialDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto nb_element = mesh.getNbElement(type, ghost_type);
const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
if (not this->rotation_matrices.exists(type, ghost_type)) {
this->rotation_matrices.alloc(0, nb_dof * nb_dof, type, ghost_type);
}
auto & rot_matrices = this->rotation_matrices(type, ghost_type);
rot_matrices.resize(nb_element);
Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
bool has_extra_normal = mesh.hasData<Real>("extra_normal", type, ghost_type);
Array<Real>::const_vector_iterator extra_normal;
if (has_extra_normal)
extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
.begin(spatial_dimension);
for (auto && tuple :
zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
make_view(rot_matrices, nb_dof, nb_dof))) {
// compute shape derivatives
auto & X = std::get<0>(tuple);
auto & R = std::get<1>(tuple);
if (has_extra_normal) {
ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
++extra_normal;
} else {
ElementClass<type>::computeRotationMatrix(
R, X, Vector<Real>(spatial_dimension));
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeShapesOnIntegrationPoints(
const Array<Real> & nodes, const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const auto & natural_coords = integration_points(type, ghost_type);
auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
auto nb_points = integration_points(type, ghost_type).cols();
auto nb_element = mesh.getNbElement(type, ghost_type);
auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
const auto dim = ElementClass<type>::getSpatialDimension();
auto itp_type = FEEngine::getInterpolationType(type);
if (not shapes.exists(itp_type, ghost_type)) {
auto size_of_shapes = this->getShapeSize(type);
this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
}
auto & shapes_ = this->shapes(itp_type, ghost_type);
shapes_.resize(nb_element * nb_points);
auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
for (auto && tuple :
zip(make_view(shapes_, nb_dof, nb_dof * nb_nodes_per_element, nb_points),
make_view(*nodes_per_element, dim, nb_nodes_per_element))) {
auto & N = std::get<0>(tuple);
auto & real_coord = std::get<1>(tuple);
ElementClass<type>::computeShapes(natural_coords, real_coord, N);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::precomputeShapeDerivativesOnIntegrationPoints(
const Array<Real> & nodes, const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const auto & natural_coords = integration_points(type, ghost_type);
const auto spatial_dimension = mesh.getSpatialDimension();
const auto natural_spatial_dimension =
ElementClass<type>::getNaturalSpaceDimension();
const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const auto nb_points = natural_coords.cols();
const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
const auto nb_element = mesh.getNbElement(type, ghost_type);
const auto nb_stress_components = ElementClass<type>::getNbStressComponents();
auto itp_type = FEEngine::getInterpolationType(type);
if (not this->shapes_derivatives.exists(itp_type, ghost_type)) {
auto size_of_shapesd = this->getShapeDerivativesSize(type);
this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
}
auto & rot_matrices = this->rotation_matrices(type, ghost_type);
Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);
auto & shapesd = this->shapes_derivatives(itp_type, ghost_type);
shapesd.resize(nb_element * nb_points);
for (auto && tuple :
zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
make_view(shapesd, nb_stress_components,
nb_nodes_per_element * nb_dof, nb_points),
make_view(rot_matrices, nb_dof, nb_dof))) {
// compute shape derivatives
auto & X = std::get<0>(tuple);
auto & B = std::get<1>(tuple);
auto & RDOFs = std::get<2>(tuple);
+
+ // /!\ This is a temporary variable with no specified size for columns !
+ // It is up to the element to resize
+ Tensor3<Real> dnds(natural_spatial_dimension,
+ ElementClass<type>::interpolation_property::dnds_columns,
+ B.size(2));
+ ElementClass<type>::computeDNDS(natural_coords, X, dnds);
+
+ Tensor3<Real> J(natural_spatial_dimension, natural_coords.rows(), natural_coords.cols());
+
+ // Computing the coordinates of the element in the natural space
auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
Matrix<Real> T(B.size(1), B.size(1));
T.block(RDOFs, 0, 0);
T.block(RDOFs, RDOFs.rows(), RDOFs.rows());
// Rotate to local basis
auto x =
(R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);
- Tensor3<Real> dnds(B.size(0), B.size(1), B.size(2));
- ElementClass<type>::computeDNDS(natural_coords, X, dnds);
-
- Tensor3<Real> J(x.rows(), natural_coords.rows(), natural_coords.cols());
-
- ElementClass<type>::computeJMat(dnds, x, J);
+ ElementClass<type>::computeJMat(natural_coords, x, J);
ElementClass<type>::computeShapeDerivatives(J, dnds, T, B);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::interpolateOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_dof,
const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof,
"The output array shape is not correct");
auto itp_type = FEEngine::getInterpolationType(type);
const auto & shapes_ = shapes(itp_type, ghost_type);
auto nb_element = mesh.getNbElement(type, ghost_type);
auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
filter_elements);
auto nb_quad_points = nb_quad_points_per_element * u_el.size();
out_uq.resize(nb_quad_points);
auto out_it = out_uq.begin_reinterpret(nb_dof, 1, nb_quad_points_per_element,
u_el.size());
auto shapes_it =
shapes_.begin_reinterpret(nb_dof, nb_dof * nb_nodes_per_element,
nb_quad_points_per_element, nb_element);
auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
nb_quad_points_per_element, u_el.size());
for_each_element(nb_element, filter_elements, [&](auto && el) {
auto & uq = *out_it;
const auto & u = *u_it;
auto N = Tensor3<Real>(shapes_it[el]);
for (auto && q : arange(uq.size(2))) {
auto uq_q = Matrix<Real>(uq(q));
auto u_q = Matrix<Real>(u(q));
auto N_q = Matrix<Real>(N(q));
uq_q.mul<false, false>(N_q, u_q);
}
++out_it;
++u_it;
});
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementKind kind>
template <ElementType type>
void ShapeStructural<kind>::gradientOnIntegrationPoints(
const Array<Real> & in_u, Array<Real> & out_nablauq, UInt nb_dof,
const GhostType & ghost_type, const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
auto itp_type = FEEngine::getInterpolationType(type);
const auto & shapesd = shapes_derivatives(itp_type, ghost_type);
auto nb_element = mesh.getNbElement(type, ghost_type);
auto element_dimension = ElementClass<type>::getSpatialDimension();
auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
filter_elements);
auto nb_quad_points = nb_quad_points_per_element * u_el.size();
out_nablauq.resize(nb_quad_points);
auto out_it = out_nablauq.begin_reinterpret(
element_dimension, 1, nb_quad_points_per_element, u_el.size());
auto shapesd_it = shapesd.begin_reinterpret(
element_dimension, nb_dof * nb_nodes_per_element,
nb_quad_points_per_element, nb_element);
auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
nb_quad_points_per_element, u_el.size());
for_each_element(nb_element, filter_elements, [&](auto && el) {
auto & nablau = *out_it;
const auto & u = *u_it;
auto B = Tensor3<Real>(shapesd_it[el]);
for (auto && q : arange(nablau.size(2))) {
auto nablau_q = Matrix<Real>(nablau(q));
auto u_q = Matrix<Real>(u(q));
auto B_q = Matrix<Real>(B(q));
nablau_q.mul<false, false>(B_q, u_q);
}
++out_it;
++u_it;
});
AKANTU_DEBUG_OUT();
}
} // namespace akantu
#endif /* __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ */
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 91f3a1581..a3f68b552 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,682 +1,683 @@
/**
* @file material.hh
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Nov 25 2015
*
* @brief Mother class for all materials
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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_factory.hh"
#include "aka_voigthelper.hh"
#include "data_accessor.hh"
#include "integration_point.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#include "mesh_events.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_HH__
#define __AKANTU_MATERIAL_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
} // namespace akantu
namespace akantu {
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element =
* ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(const ElementType & el_type,
* Array<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : public Memory,
public DataAccessor<Element>,
public Parsable,
public MeshEventHandler,
protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
#if __cplusplus > 199711L
Material(const Material & mat) = delete;
Material & operator=(const Material & mat) = delete;
#endif
/// Initialize material with defaults
Material(SolidMechanicsModel & model, const ID & id = "");
/// Initialize material with custom mesh & fe_engine
Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id = "");
/// Destructor
~Material() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(__attribute__((unused))
const ElementType & el_type,
__attribute__((unused))
Array<Real> & tangent_matrix,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type,
GhostType ghost_type = _not_ghost);
/// compute the potential energy for an element
virtual void
computePotentialEnergyByElement(__attribute__((unused)) ElementType type,
__attribute__((unused)) UInt index,
__attribute__((unused))
Vector<Real> & epot_on_quad_points) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
virtual void updateEnergies(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
virtual void updateEnergiesAfterDamage(__attribute__((unused))
ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
/// set the material to steady state (to be implemented for materials that
/// need it)
virtual void setToSteadyState(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
/// function called to update the internal parameters when the modifiable
/// parameters are modified
virtual void updateInternalParameters() {}
public:
/// extrapolate internal values
virtual void extrapolateInternal(const ID & id, const Element & element,
const Matrix<Real> & points,
Matrix<Real> & extrapolated);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(__attribute__((unused))
const Element & element) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(__attribute__((unused))
const Element & element) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// get a material celerity to compute the stable time step (default: is the
/// push wave speed)
virtual Real getCelerity(const Element & element) const {
return getPushWaveSpeed(element);
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <typename T>
void registerInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
template <typename T>
void unregisterInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
// virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleInternalForces(GhostType ghost_type);
/// save the stress in the previous_stress if needed
virtual void savePreviousState();
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
// virtual void
// computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// add an element to the local mesh filter
inline UInt addElement(const ElementType & type, UInt element,
const GhostType & ghost_type);
inline UInt addElement(const Element & element);
/// add many elements at once
void addElements(const Array<Element> & elements_to_add);
/// remove many element at once
void removeElements(const Array<Element> & elements_to_remove);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray<Real> & result,
const GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
ElementTypeMapArray<Real> & by_elem_result,
const GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(
const ElementTypeMapArray<Real> & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/* ------------------------------------------------------------------------ */
inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const;
/// compute the potential energy by element
void computePotentialEnergyByElements();
/// resize the intenals arrays
virtual void resizeInternals();
/* ------------------------------------------------------------------------ */
/* Finite deformation functions */
/* This functions area implementing what is described in the paper of Bathe */
/* et al, in IJNME, Finite Element Formulations For Large Deformation */
/* Dynamic Analysis, Vol 9, 353-386, 1975 */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the residual
template <UInt dim> void assembleInternalForces(GhostType ghost_type);
/// Computation of Cauchy stress tensor in the case of finite deformation from
/// the 2nd Piola-Kirchhoff for a given element type
template <UInt dim>
void computeCauchyStress(ElementType el_type,
GhostType ghost_type = _not_ghost);
/// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
/// gradient
template <UInt dim>
inline void computeCauchyStressOnQuad(const Matrix<Real> & F,
const Matrix<Real> & S,
Matrix<Real> & cauchy,
const Real & C33 = 1.0) const;
template <UInt dim>
void computeAllStressesFromTangentModuli(const ElementType & type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type);
/// assembling in finite deformation
template <UInt dim>
void assembleStiffnessMatrixNL(const ElementType & type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrixL2(const ElementType & type,
GhostType ghost_type);
/// Size of the Stress matrix for the case of finite deformation see: Bathe et
/// al, IJNME, Vol 9, 353-386, 1975
inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const;
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
/// 1975
template <UInt dim>
inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & sigma);
/// write the stress tensor in the Voigt notation.
template <UInt dim>
inline void setCauchyStressArray(const Matrix<Real> & S_t,
Matrix<Real> & sigma_voight);
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
public:
template <UInt dim>
static inline void gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F);
static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
static inline void leftCauchy(const Matrix<Real> & F, Matrix<Real> & B);
template <UInt dim>
static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline void gradUToGreenStrain(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
static inline Real stressToVonMises(const Matrix<Real> & stress);
protected:
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/// converts global quadrature point to local quadrature point
inline IntegrationPoint
convertToLocalPoint(const IntegrationPoint & global_point) const;
/// converts local quadrature point to global quadrature point
inline IntegrationPoint
convertToGlobalPoint(const IntegrationPoint & local_point) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) override;
template <typename T>
inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID()) const;
template <typename T>
inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
void onNodesAdded(const Array<UInt> &, const NewNodesEvent &) override{};
void onNodesRemoved(const Array<UInt> &, const Array<UInt> &,
const RemovedNodesEvent &) override{};
void onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) override;
void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) override;
void onElementsChanged(const Array<Element> &, const Array<Element> &,
const ElementTypeMapArray<UInt> &,
const ChangedElementsEvent &) override{};
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void beforeSolveStep();
virtual void afterSolveStep();
void onDamageIteration() override;
void onDamageUpdate() override;
void onDump() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Name, name, const std::string &);
AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(ID, Memory::getID(), const ID &);
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// return the potential energy for the subset of elements contained by the
/// material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(ElementType & type, UInt index);
/// return the energy (identified by id) for the subset of elements contained
/// by the material
virtual Real getEnergy(const std::string & energy_id);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(const std::string & energy_id, ElementType type,
UInt index);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(ElementFilter, element_filter,
const ElementTypeMapArray<UInt> &);
AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
bool isNonLocal() const { return is_non_local; }
template <typename T>
const Array<T> & getArray(const ID & id, const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
template <typename T>
Array<T> & getArray(const ID & id, const ElementType & type,
const GhostType & ghost_type = _not_ghost);
template <typename T>
const InternalField<T> & getInternal(const ID & id) const;
template <typename T> InternalField<T> & getInternal(const ID & id);
template <typename T>
inline bool isInternal(const ID & id, const ElementKind & element_kind) const;
template <typename T>
ElementTypeMap<UInt>
getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const;
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
template <typename T> inline void setParam(const ID & param, T value);
+ template <typename T> inline void setParamNoUpdate(const ID & param, T value);
inline const Parameter & getParam(const ID & param) const;
template <typename T>
void flattenInternal(const std::string & field_id,
ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) const;
/// apply a constant eigengrad_u everywhere in the material
virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
const GhostType = _not_ghost);
/// specify if the matrix need to be recomputed for this material
virtual bool hasStiffnessMatrixChanged() { return true; }
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// boolean to know if the material has been initialized
bool is_init;
std::map<ID, InternalField<Real> *> internal_vectors_real;
std::map<ID, InternalField<UInt> *> internal_vectors_uint;
std::map<ID, InternalField<bool> *> internal_vectors_bool;
protected:
/// Link to the fem object in the model
FEEngine & fem;
/// Finite deformation
bool finite_deformation;
/// Finite deformation
bool inelastic_deformation;
/// material name
std::string name;
/// The model to witch the material belong
SolidMechanicsModel & model;
/// density : rho
Real rho;
/// spatial dimension
UInt spatial_dimension;
/// list of element handled by the material
ElementTypeMapArray<UInt> element_filter;
/// stresses arrays ordered by element types
InternalField<Real> stress;
/// eigengrad_u arrays ordered by element types
InternalField<Real> eigengradu;
/// grad_u arrays ordered by element types
InternalField<Real> gradu;
/// Green Lagrange strain (Finite deformation)
InternalField<Real> green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types
/// (Finite deformation)
InternalField<Real> piola_kirchhoff_2;
/// potential energy by element
InternalField<Real> potential_energy;
/// tell if using in non local mode or not
bool is_non_local;
/// tell if the material need the previous stress state
bool use_previous_stress;
/// tell if the material need the previous strain state
bool use_previous_gradu;
/// elemental field interpolation coordinates
InternalField<Real> interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField<Real> interpolation_points_matrices;
/// vector that contains the names of all the internals that need to
/// be transferred when material interfaces move
std::vector<ID> internals_to_transfer;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Material & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "material_inline_impl.cc"
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeStress. This macro in addition to write the loop
/// provides two tensors (matrices) sigma and grad_u
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type) \
.end(this->spatial_dimension, this->spatial_dimension); \
\
this->stress(el_type, ghost_type) \
.resize(this->gradu(el_type, ghost_type).size()); \
\
Array<Real>::iterator<Matrix<Real>> stress_it = \
this->stress(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
\
if (this->isFiniteDeformation()) { \
this->piola_kirchhoff_2(el_type, ghost_type) \
.resize(this->gradu(el_type, ghost_type).size()); \
stress_it = this->piola_kirchhoff_2(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
} \
\
for (; gradu_it != gradu_end; ++gradu_it, ++stress_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma = *stress_it
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeTangentModuli. This macro in addition to write the
/// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
/// where the elemental tangent moduli should be stored in Voigt Notation
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type) \
.end(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator sigma_it = \
this->stress(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
\
tangent_mat.resize(this->gradu(el_type, ghost_type).size()); \
\
UInt tangent_size = \
this->getTangentStiffnessVoigtSize(this->spatial_dimension); \
Array<Real>::matrix_iterator tangent_it = \
tangent_mat.begin(tangent_size, tangent_size); \
\
for (; gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma_tensor = *sigma_it; \
Matrix<Real> & tangent = *tangent_it
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
/* -------------------------------------------------------------------------- */
namespace akantu {
using MaterialFactory =
Factory<Material, ID, UInt, const ID &, SolidMechanicsModel &, const ID &>;
} // namespace akantu
#define INSTANTIATE_MATERIAL_ONLY(mat_name) \
template class mat_name<1>; \
template class mat_name<2>; \
template class mat_name<3>
#define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \
[](UInt dim, const ID &, SolidMechanicsModel & model, \
const ID & id) -> std::unique_ptr<Material> { \
switch (dim) { \
case 1: \
return std::make_unique<mat_name<1>>(model, id); \
case 2: \
return std::make_unique<mat_name<2>>(model, id); \
case 3: \
return std::make_unique<mat_name<3>>(model, id); \
default: \
AKANTU_EXCEPTION("The dimension " \
<< dim << "is not a valid dimension for the material " \
<< #id); \
} \
}
#define INSTANTIATE_MATERIAL(id, mat_name) \
INSTANTIATE_MATERIAL_ONLY(mat_name); \
static bool material_is_alocated_##id [[gnu::unused]] = \
MaterialFactory::getInstance().registerAllocator( \
#id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name))
#endif /* __AKANTU_MATERIAL_HH__ */
diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc
index 4b06b4d13..1f6d1c0dd 100644
--- a/src/model/solid_mechanics/material_inline_impl.cc
+++ b/src/model/solid_mechanics/material_inline_impl.cc
@@ -1,464 +1,475 @@
/**
* @file material_inline_impl.cc
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Wed Nov 25 2015
*
* @brief Implementation of the inline functions of the class material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_INLINE_IMPL_CC__
#define __AKANTU_MATERIAL_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
inline UInt Material::addElement(const ElementType & type, UInt element,
const GhostType & ghost_type) {
Array<UInt> & el_filter = this->element_filter(type, ghost_type);
el_filter.push_back(element);
return el_filter.size() - 1;
}
/* -------------------------------------------------------------------------- */
inline UInt Material::addElement(const Element & element) {
return this->addElement(element.type, element.element, element.ghost_type);
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const {
return (dim * (dim - 1) / 2 + dim);
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getCauchyStressMatrixSize(UInt dim) const {
return (dim * dim);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline void Material::gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F) {
AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim * dim,
"The dimension of the tensor F should be greater or "
"equal to the dimension of the tensor grad_u.");
F.eye();
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
F(i, j) += grad_u(i, j);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline void Material::computeCauchyStressOnQuad(const Matrix<Real> & F,
const Matrix<Real> & piola,
Matrix<Real> & sigma,
const Real & C33) const {
Real J = F.det() * sqrt(C33);
Matrix<Real> F_S(dim, dim);
F_S.mul<false, false>(F, piola);
Real constant = J ? 1. / J : 0;
sigma.mul<false, true>(F_S, F, constant);
}
/* -------------------------------------------------------------------------- */
inline void Material::rightCauchy(const Matrix<Real> & F, Matrix<Real> & C) {
C.mul<true, false>(F, F);
}
/* -------------------------------------------------------------------------- */
inline void Material::leftCauchy(const Matrix<Real> & F, Matrix<Real> & B) {
B.mul<false, true>(F, F);
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline void Material::gradUToEpsilon(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon) {
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i));
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline void Material::gradUToGreenStrain(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon) {
epsilon.mul<true, false>(grad_u, grad_u, .5);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) += 0.5 * (grad_u(i, j) + grad_u(j, i));
}
/* -------------------------------------------------------------------------- */
inline Real Material::stressToVonMises(const Matrix<Real> & stress) {
// compute deviatoric stress
UInt dim = stress.cols();
Matrix<Real> deviatoric_stress =
Matrix<Real>::eye(dim, -1. * stress.trace() / 3.);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
deviatoric_stress(i, j) += stress(i, j);
// return Von Mises stress
return std::sqrt(3. * deviatoric_stress.doubleDot(deviatoric_stress) / 2.);
}
/* ---------------------------------------------------------------------------*/
template <UInt dim>
inline void Material::setCauchyStressArray(const Matrix<Real> & S_t,
Matrix<Real> & sigma_voight) {
AKANTU_DEBUG_IN();
sigma_voight.clear();
// see Finite element formulations for large deformation dynamic analysis,
// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau
/*
* 1d: [ s11 ]'
* 2d: [ s11 s22 s12 ]'
* 3d: [ s11 s22 s33 s23 s13 s12 ]
*/
for (UInt i = 0; i < dim; ++i) // diagonal terms
sigma_voight(i, 0) = S_t(i, i);
for (UInt i = 1; i < dim; ++i) // term s12 in 2D and terms s23 s13 in 3D
sigma_voight(dim + i - 1, 0) = S_t(dim - i - 1, dim - 1);
for (UInt i = 2; i < dim; ++i) // term s13 in 3D
sigma_voight(dim + i, 0) = S_t(0, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
inline void Material::setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & sigma) {
AKANTU_DEBUG_IN();
sigma.clear();
/// see Finite ekement formulations for large deformation dynamic analysis,
/// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau
for (UInt i = 0; i < dim; ++i) {
for (UInt m = 0; m < dim; ++m) {
for (UInt n = 0; n < dim; ++n) {
sigma(i * dim + m, i * dim + n) = S_t(m, n);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline Element
Material::convertToLocalElement(const Element & global_element) const {
UInt ge = global_element.element;
#ifndef AKANTU_NDEBUG
UInt model_mat_index = this->model.getMaterialByElement(
global_element.type, global_element.ghost_type)(ge);
UInt mat_index = this->model.getMaterialIndex(this->name);
AKANTU_DEBUG_ASSERT(model_mat_index == mat_index,
"Conversion of a global element in a local element for "
"the wrong material "
<< this->name << std::endl);
#endif
UInt le = this->model.getMaterialLocalNumbering(
global_element.type, global_element.ghost_type)(ge);
Element tmp_quad{global_element.type, le, global_element.ghost_type};
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline Element
Material::convertToGlobalElement(const Element & local_element) const {
UInt le = local_element.element;
UInt ge =
this->element_filter(local_element.type, local_element.ghost_type)(le);
Element tmp_quad{local_element.type, ge, local_element.ghost_type};
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline IntegrationPoint
Material::convertToLocalPoint(const IntegrationPoint & global_point) const {
const FEEngine & fem = this->model.getFEEngine();
UInt nb_quad = fem.getNbIntegrationPoints(global_point.type);
Element el =
this->convertToLocalElement(static_cast<const Element &>(global_point));
IntegrationPoint tmp_quad(el, global_point.num_point, nb_quad);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline IntegrationPoint
Material::convertToGlobalPoint(const IntegrationPoint & local_point) const {
const FEEngine & fem = this->model.getFEEngine();
UInt nb_quad = fem.getNbIntegrationPoints(local_point.type);
Element el =
this->convertToGlobalElement(static_cast<const Element &>(local_point));
IntegrationPoint tmp_quad(el, local_point.num_point, nb_quad);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const {
if (tag == _gst_smm_stress) {
return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension *
spatial_dimension * sizeof(Real) *
this->getModel().getNbIntegrationPoints(elements);
}
return 0;
}
/* -------------------------------------------------------------------------- */
inline void Material::packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const {
if (tag == _gst_smm_stress) {
if (this->isFiniteDeformation()) {
packElementDataHelper(piola_kirchhoff_2, buffer, elements);
packElementDataHelper(gradu, buffer, elements);
}
packElementDataHelper(stress, buffer, elements);
}
}
/* -------------------------------------------------------------------------- */
inline void Material::unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) {
if (tag == _gst_smm_stress) {
if (this->isFiniteDeformation()) {
unpackElementDataHelper(piola_kirchhoff_2, buffer, elements);
unpackElementDataHelper(gradu, buffer, elements);
}
unpackElementDataHelper(stress, buffer, elements);
}
}
/* -------------------------------------------------------------------------- */
inline const Parameter & Material::getParam(const ID & param) const {
try {
return get(param);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the material "
<< getID());
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Material::setParam(const ID & param, T value) {
try {
set<T>(param, value);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the material "
<< getID());
}
updateInternalParameters();
}
+/* -------------------------------------------------------------------------- */
+template <typename T>
+inline void Material::setParamNoUpdate(const ID & param, T value) {
+ try {
+ set<T>(param, value);
+ } catch (...) {
+ AKANTU_EXCEPTION("No parameter " << param << " in the material "
+ << getID());
+ }
+}
+
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Material::packElementDataHelper(
const ElementTypeMapArray<T> & data_to_pack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) const {
DataAccessor::packElementalDataHelper<T>(data_to_pack, buffer, elements, true,
model.getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Material::unpackElementDataHelper(
ElementTypeMapArray<T> & data_to_unpack, CommunicationBuffer & buffer,
const Array<Element> & elements, const ID & fem_id) {
DataAccessor::unpackElementalDataHelper<T>(data_to_unpack, buffer, elements,
true, model.getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template <>
inline void Material::registerInternal<Real>(InternalField<Real> & vect) {
internal_vectors_real[vect.getID()] = &vect;
}
template <>
inline void Material::registerInternal<UInt>(InternalField<UInt> & vect) {
internal_vectors_uint[vect.getID()] = &vect;
}
template <>
inline void Material::registerInternal<bool>(InternalField<bool> & vect) {
internal_vectors_bool[vect.getID()] = &vect;
}
/* -------------------------------------------------------------------------- */
template <>
inline void Material::unregisterInternal<Real>(InternalField<Real> & vect) {
internal_vectors_real.erase(vect.getID());
}
template <>
inline void Material::unregisterInternal<UInt>(InternalField<UInt> & vect) {
internal_vectors_uint.erase(vect.getID());
}
template <>
inline void Material::unregisterInternal<bool>(InternalField<bool> & vect) {
internal_vectors_bool.erase(vect.getID());
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline bool Material::isInternal(__attribute__((unused)) const ID & id,
__attribute__((unused))
const ElementKind & element_kind) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
template <>
inline bool Material::isInternal<Real>(const ID & id,
const ElementKind & element_kind) const {
auto internal_array = internal_vectors_real.find(this->getID() + ":" + id);
if (internal_array == internal_vectors_real.end() ||
internal_array->second->getElementKind() != element_kind)
return false;
return true;
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline ElementTypeMap<UInt>
Material::getInternalDataPerElem(const ID & field_id,
const ElementKind & element_kind) const {
if (!this->template isInternal<T>(field_id, element_kind))
AKANTU_EXCEPTION("Cannot find internal field " << id << " in material "
<< this->name);
const InternalField<T> & internal_field =
this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
UInt nb_data_per_quad = internal_field.getNbComponent();
ElementTypeMap<UInt> res;
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
using type_iterator = typename InternalField<T>::type_iterator;
type_iterator tit = internal_field.firstType(*gt);
type_iterator tend = internal_field.lastType(*gt);
for (; tit != tend; ++tit) {
UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(*tit, *gt);
res(*tit, *gt) = nb_data_per_quad * nb_quadrature_points;
}
}
return res;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Material::flattenInternal(const std::string & field_id,
ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type,
ElementKind element_kind) const {
if (!this->template isInternal<T>(field_id, element_kind))
AKANTU_EXCEPTION("Cannot find internal field " << id << " in material "
<< this->name);
const InternalField<T> & internal_field =
this->template getInternal<T>(field_id);
const FEEngine & fe_engine = internal_field.getFEEngine();
const Mesh & mesh = fe_engine.getMesh();
using type_iterator = typename InternalField<T>::filter_type_iterator;
type_iterator tit = internal_field.filterFirstType(ghost_type);
type_iterator tend = internal_field.filterLastType(ghost_type);
for (; tit != tend; ++tit) {
ElementType type = *tit;
const Array<Real> & src_vect = internal_field(type, ghost_type);
const Array<UInt> & filter = internal_field.getFilter(type, ghost_type);
// total number of elements in the corresponding mesh
UInt nb_element_dst = mesh.getNbElement(type, ghost_type);
// number of element in the internal field
UInt nb_element_src = filter.size();
// number of quadrature points per elem
UInt nb_quad_per_elem = fe_engine.getNbIntegrationPoints(type);
// number of data per quadrature point
UInt nb_data_per_quad = internal_field.getNbComponent();
if (!internal_flat.exists(type, ghost_type)) {
internal_flat.alloc(nb_element_dst * nb_quad_per_elem, nb_data_per_quad,
type, ghost_type);
}
if (nb_element_src == 0)
continue;
// number of data per element
UInt nb_data = nb_quad_per_elem * nb_data_per_quad;
Array<Real> & dst_vect = internal_flat(type, ghost_type);
dst_vect.resize(nb_element_dst * nb_quad_per_elem);
Array<UInt>::const_scalar_iterator it = filter.begin();
Array<UInt>::const_scalar_iterator end = filter.end();
Array<Real>::const_vector_iterator it_src =
src_vect.begin_reinterpret(nb_data, nb_element_src);
Array<Real>::vector_iterator it_dst =
dst_vect.begin_reinterpret(nb_data, nb_element_dst);
for (; it != end; ++it, ++it_src) {
it_dst[*it] = *it_src;
}
}
}
} // akantu
#endif /* __AKANTU_MATERIAL_INLINE_IMPL_CC__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
index e24a266cd..5f08a8369 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
@@ -1,308 +1,266 @@
/**
* @file material_elastic_linear_anisotropic.cc
*
* @author Till Junge <till.junge@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Sep 25 2013
* @date last modification: Thu Oct 15 2015
*
* @brief Anisotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2014, 2015 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 "material_elastic_linear_anisotropic.hh"
#include "solid_mechanics_model.hh"
#include <algorithm>
#include <sstream>
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt dim>
MaterialElasticLinearAnisotropic<dim>::MaterialElasticLinearAnisotropic(
SolidMechanicsModel & model, const ID & id, bool symmetric)
: Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim),
C(voigt_h::size, voigt_h::size), eigC(voigt_h::size),
symmetric(symmetric), alpha(0), was_stiffness_assembled(false) {
AKANTU_DEBUG_IN();
this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
(*this->dir_vecs.back())[0] = 1.;
this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of main material axis");
if (dim > 1) {
this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
(*this->dir_vecs.back())[1] = 1.;
this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of secondary material axis");
}
if (dim > 2) {
this->dir_vecs.push_back(std::make_unique<Vector<Real>>(dim));
(*this->dir_vecs.back())[2] = 1.;
this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of tertiary material axis");
}
for (UInt i = 0; i < voigt_h::size; ++i) {
UInt start = 0;
if (this->symmetric) {
start = i;
}
for (UInt j = start; j < voigt_h::size; ++j) {
std::stringstream param("C");
param << "C" << i + 1 << j + 1;
this->registerParam(param.str(), this->Cprime(i, j), Real(0.),
_pat_parsmod, "Coefficient " + param.str());
}
}
this->registerParam("alpha", this->alpha, _pat_parsmod,
"Proportion of viscous stress");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim> void MaterialElasticLinearAnisotropic<dim>::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
updateInternalParameters();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
void MaterialElasticLinearAnisotropic<dim>::updateInternalParameters() {
Material::updateInternalParameters();
if (this->symmetric) {
for (UInt i = 0; i < voigt_h::size; ++i) {
for (UInt j = i + 1; j < voigt_h::size; ++j) {
this->Cprime(j, i) = this->Cprime(i, j);
}
}
}
this->rotateCprime();
this->C.eig(this->eigC);
this->was_stiffness_assembled = false;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void MaterialElasticLinearAnisotropic<Dim>::rotateCprime() {
// start by filling the empty parts fo Cprime
UInt diff = Dim * Dim - voigt_h::size;
for (UInt i = voigt_h::size; i < Dim * Dim; ++i) {
for (UInt j = 0; j < Dim * Dim; ++j) {
this->Cprime(i, j) = this->Cprime(i - diff, j);
}
}
for (UInt i = 0; i < Dim * Dim; ++i) {
for (UInt j = voigt_h::size; j < Dim * Dim; ++j) {
this->Cprime(i, j) = this->Cprime(i, j - diff);
}
}
// construction of rotator tensor
// normalise rotation matrix
for (UInt j = 0; j < Dim; ++j) {
Vector<Real> rot_vec = this->rot_mat(j);
rot_vec = *this->dir_vecs[j];
rot_vec.normalize();
}
// make sure the vectors form a right-handed base
Vector<Real> test_axis(3);
Vector<Real> v1(3), v2(3), v3(3);
if (Dim == 2) {
for (UInt i = 0; i < Dim; ++i) {
v1[i] = this->rot_mat(0, i);
v2[i] = this->rot_mat(1, i);
v3[i] = 0.;
}
v3[2] = 1.;
v1[2] = 0.;
v2[2] = 0.;
} else if (Dim == 3) {
v1 = this->rot_mat(0);
v2 = this->rot_mat(1);
v3 = this->rot_mat(2);
}
test_axis.crossProduct(v1, v2);
test_axis -= v3;
if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
AKANTU_DEBUG_ERROR("The axis vectors do not form a right-handed coordinate "
<< "system. I. e., ||n1 x n2 - n3|| should be zero, but "
<< "it is " << test_axis.norm() << ".");
}
// create the rotator and the reverse rotator
Matrix<Real> rotator(Dim * Dim, Dim * Dim);
Matrix<Real> revrotor(Dim * Dim, Dim * Dim);
for (UInt i = 0; i < Dim; ++i) {
for (UInt j = 0; j < Dim; ++j) {
for (UInt k = 0; k < Dim; ++k) {
for (UInt l = 0; l < Dim; ++l) {
UInt I = voigt_h::mat[i][j];
UInt J = voigt_h::mat[k][l];
rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
}
}
}
}
// create the full rotated matrix
Matrix<Real> Cfull(Dim * Dim, Dim * Dim);
Cfull = rotator * Cprime * revrotor;
for (UInt i = 0; i < voigt_h::size; ++i) {
for (UInt j = 0; j < voigt_h::size; ++j) {
this->C(i, j) = Cfull(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
-void MaterialElasticLinearAnisotropic<dim>::computeStress(
- ElementType el_type, GhostType ghost_type) {
+void MaterialElasticLinearAnisotropic<dim>::computeStress(ElementType el_type,
+ GhostType ghost_type) {
// Wikipedia convention:
// 2*eps_ij (i!=j) = voigt_eps_I
// http://en.wikipedia.org/wiki/Voigt_notation
AKANTU_DEBUG_IN();
- Array<Real>::iterator<Matrix<Real>> gradu_it =
- this->gradu(el_type, ghost_type).begin(dim, dim);
- Array<Real>::iterator<Matrix<Real>> gradu_end =
- this->gradu(el_type, ghost_type).end(dim, dim);
-
- UInt nb_quad_pts = gradu_end - gradu_it;
-
- // create array for strains and stresses of all dof of all gauss points
- // for efficient computation of stress
- Matrix<Real> voigt_strains(voigt_h::size, nb_quad_pts);
- Matrix<Real> voigt_stresses(voigt_h::size, nb_quad_pts);
-
- // copy strains
Matrix<Real> strain(dim, dim);
- for (UInt q = 0; gradu_it != gradu_end; ++gradu_it, ++q) {
- Matrix<Real> & grad_u = *gradu_it;
-
- for (UInt I = 0; I < voigt_h::size; ++I) {
- Real voigt_factor = voigt_h::factors[I];
- UInt i = voigt_h::vec[I][0];
- UInt j = voigt_h::vec[I][1];
-
- voigt_strains(I, q) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.;
- }
- }
-
- // compute the strain rate proportional part if needed
- // bool viscous = this->alpha == 0.; // only works if default value
- bool viscous = false;
- if (viscous) {
- Array<Real> strain_rate(0, dim * dim, "strain_rate");
-
- Array<Real> & velocity = this->model.getVelocity();
- const Array<UInt> & elem_filter = this->element_filter(el_type, ghost_type);
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
- this->fem.gradientOnIntegrationPoints(velocity, strain_rate, dim, el_type,
- ghost_type, elem_filter);
+ this->computeStressOnQuad(strain,sigma);
- Array<Real>::matrix_iterator gradu_dot_it = strain_rate.begin(dim, dim);
- Array<Real>::matrix_iterator gradu_dot_end = strain_rate.end(dim, dim);
-
- Matrix<Real> strain_dot(dim, dim);
- for (UInt q = 0; gradu_dot_it != gradu_dot_end; ++gradu_dot_it, ++q) {
- Matrix<Real> & grad_u_dot = *gradu_dot_it;
-
- for (UInt I = 0; I < voigt_h::size; ++I) {
- Real voigt_factor = voigt_h::factors[I];
- UInt i = voigt_h::vec[I][0];
- UInt j = voigt_h::vec[I][1];
-
- voigt_strains(I, q) = this->alpha * voigt_factor *
- (grad_u_dot(i, j) + grad_u_dot(j, i)) / 2.;
- }
- }
- }
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
+}
- // compute stresses
- voigt_stresses = this->C * voigt_strains;
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+void MaterialElasticLinearAnisotropic<dim>::computeTangentModuli(
+ const ElementType & el_type, Array<Real> & tangent_matrix,
+ GhostType ghost_type) {
+ AKANTU_DEBUG_IN();
- // copy stresses back
- Array<Real>::iterator<Matrix<Real>> stress_it =
- this->stress(el_type, ghost_type).begin(dim, dim);
+ MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
- Array<Real>::iterator<Matrix<Real>> stress_end =
- this->stress(el_type, ghost_type).end(dim, dim);
+ this->computeTangentModuliOnQuad(tangent);
- for (UInt q = 0; stress_it != stress_end; ++stress_it, ++q) {
- Matrix<Real> & stress = *stress_it;
+ MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
- for (UInt I = 0; I < voigt_h::size; ++I) {
- UInt i = voigt_h::vec[I][0];
- UInt j = voigt_h::vec[I][1];
- stress(i, j) = stress(j, i) = voigt_stresses(I, q);
- }
- }
+ this->was_stiffness_assembled = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
-void MaterialElasticLinearAnisotropic<dim>::computeTangentModuli(
- const ElementType & el_type, Array<Real> & tangent_matrix,
- GhostType ghost_type) {
+void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergy(
+ ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
- MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
- tangent.copy(this->C);
- MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
+ Material::computePotentialEnergy(el_type, ghost_type);
- this->was_stiffness_assembled = true;
+ AKANTU_DEBUG_ASSERT(!this->finite_deformation,
+ "finite deformation not possible in material anisotropic "
+ "(TO BE IMPLEMENTED)");
+
+ if (ghost_type != _not_ghost)
+ return;
+ Array<Real>::scalar_iterator epot =
+ this->potential_energy(el_type, ghost_type).begin();
+
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
+
+ computePotentialEnergyOnQuad(grad_u, sigma, *epot);
+ ++epot;
+
+ MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt dim>
Real MaterialElasticLinearAnisotropic<dim>::getCelerity(
__attribute__((unused)) const Element & element) const {
return std::sqrt(this->eigC(0) / rho);
}
/* -------------------------------------------------------------------------- */
INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic);
} // akantu
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
index dec6765fd..2cfa119bb 100644
--- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
@@ -1,122 +1,139 @@
/**
* @file material_elastic_linear_anisotropic.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Tue Aug 18 2015
*
* @brief Orthotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "material.hh"
#include "material_elastic.hh"
/* -------------------------------------------------------------------------- */
#include <vector>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH__
#define __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH__
namespace akantu {
/**
* General linear anisotropic elastic material
* The only constraint on the elastic tensor is that it can be represented
* as a symmetric 6x6 matrix (3D) or 3x3 (2D).
*
* parameters in the material files :
* - rho : density (default: 0)
* - C_ij : entry on the stiffness
*/
template <UInt Dim> class MaterialElasticLinearAnisotropic : public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElasticLinearAnisotropic(SolidMechanicsModel & model,
const ID & id = "", bool symmetric = true);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial() override;
/// constitutive law for all element of a type
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override;
/// compute the tangent stiffness matrix for an element type
void computeTangentModuli(const ElementType & el_type,
Array<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost) override;
+ /// compute the elastic potential energy
+ void computePotentialEnergy(ElementType el_type,
+ GhostType ghost_type = _not_ghost) override;
+
void updateInternalParameters() override;
bool hasStiffnessMatrixChanged() override {
return (!was_stiffness_assembled);
}
protected:
// compute C from Cprime
void rotateCprime();
+ /// constitutive law for a given quadrature point
+ inline void computeStressOnQuad(const Matrix<Real> & grad_u,
+ Matrix<Real> & sigma) const;
+
+ /// tangent matrix for a given quadrature point
+ inline void computeTangentModuliOnQuad(Matrix<Real> & tangent) const;
+
+ inline void computePotentialEnergyOnQuad(const Matrix<Real> & grad_u,
+ const Matrix<Real> & sigma,
+ Real & epot);
+
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// compute max wave celerity
Real getCelerity(const Element & element) const override;
AKANTU_GET_MACRO(VoigtStiffness, C, Matrix<Real>);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
using voigt_h = VoigtHelper<Dim>;
/// direction matrix and vectors
std::vector<std::unique_ptr<Vector<Real>>> dir_vecs;
Matrix<Real> rot_mat;
/// Elastic stiffness tensor in material frame and full vectorised notation
Matrix<Real> Cprime;
/// Elastic stiffness tensor in voigt notation
Matrix<Real> C;
/// eigenvalues of stiffness tensor
Vector<Real> eigC;
bool symmetric;
/// viscous proportion
Real alpha;
/// defines if the stiffness was computed
bool was_stiffness_assembled;
};
} // akantu
+#include "material_elastic_linear_anisotropic_inline_impl.cc"
+
#endif /* __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.cc
new file mode 100644
index 000000000..900b8d5ae
--- /dev/null
+++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.cc
@@ -0,0 +1,96 @@
+/**
+ * @file material_elastic_linear_anisotropic_inline_impl.cc
+ *
+ * @author Enrico Milanese <enrico.milanese@epfl.ch>
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @date creation: Mon Feb 12 2018
+ * @date last modification: Mon Feb 12 2018
+ *
+ * @brief Implementation of the inline functions of the material elastic linear
+ * anisotropic
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2012, 2014, 2015 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 "material_elastic_linear_anisotropic.hh"
+/* -------------------------------------------------------------------------- */
+
+#ifndef __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_INLINE_IMPL_CC__
+#define __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_INLINE_IMPL_CC__
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+inline void MaterialElasticLinearAnisotropic<dim>::computeStressOnQuad(const Matrix<Real> & grad_u,
+ Matrix<Real> & sigma) const {
+ // Wikipedia convention:
+ // 2*eps_ij (i!=j) = voigt_eps_I
+ // http://en.wikipedia.org/wiki/Voigt_notation
+ Vector<Real> voigt_strain(voigt_h::size);
+ Vector<Real> voigt_stress(voigt_h::size);
+
+ for (UInt I = 0; I < voigt_h::size; ++I){
+ Real voigt_factor = voigt_h::factors[I];
+ UInt i = voigt_h::vec[I][0];
+ UInt j = voigt_h::vec[I][1];
+
+ voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.;
+
+ }
+
+ voigt_stress = this->C * voigt_strain;
+
+ for (UInt I = 0; I < voigt_h::size; ++I){
+ UInt i = voigt_h::vec[I][0];
+ UInt j = voigt_h::vec[I][1];
+
+ sigma(i, j) = sigma(j, i) = voigt_stress(I);
+
+ }
+
+}
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+inline void MaterialElasticLinearAnisotropic<dim>::computeTangentModuliOnQuad(Matrix<Real> & tangent) const {
+
+ tangent.copy(this->C);
+
+}
+
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+inline void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergyOnQuad(
+ const Matrix<Real> & grad_u, const Matrix<Real> & sigma, Real & epot) {
+
+ AKANTU_DEBUG_ASSERT(this->symmetric,
+ "The elastic constants matrix is not symmetric,"
+ "energy is not path independent.");
+
+ epot = .5 * sigma.doubleDot(grad_u);
+}
+
+} // akantu
+
+#endif /* __AKANTU_MATERIAL_ELASTIC_LINEAR_ANISOTROPIC_INLINE_IMPL_CC__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc
index 22815792a..331a25fbe 100644
--- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc
@@ -1,212 +1,168 @@
/**
* @file material_elastic_orthotropic.cc
*
* @author Till Junge <till.junge@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Thu Oct 15 2015
*
* @brief Orthotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "material_elastic_orthotropic.hh"
#include "solid_mechanics_model.hh"
#include <algorithm>
namespace akantu {
/* -------------------------------------------------------------------------- */
template <UInt Dim>
MaterialElasticOrthotropic<Dim>::MaterialElasticOrthotropic(
SolidMechanicsModel & model, const ID & id)
: MaterialElasticLinearAnisotropic<Dim>(model, id) {
AKANTU_DEBUG_IN();
this->registerParam("E1", E1, Real(0.), _pat_parsmod, "Young's modulus (n1)");
this->registerParam("E2", E2, Real(0.), _pat_parsmod, "Young's modulus (n2)");
this->registerParam("nu12", nu12, Real(0.), _pat_parsmod,
"Poisson's ratio (12)");
this->registerParam("G12", G12, Real(0.), _pat_parsmod, "Shear modulus (12)");
if (Dim > 2) {
this->registerParam("E3", E3, Real(0.), _pat_parsmod,
"Young's modulus (n3)");
this->registerParam("nu13", nu13, Real(0.), _pat_parsmod,
"Poisson's ratio (13)");
this->registerParam("nu23", nu23, Real(0.), _pat_parsmod,
"Poisson's ratio (23)");
this->registerParam("G13", G13, Real(0.), _pat_parsmod,
"Shear modulus (13)");
this->registerParam("G23", G23, Real(0.), _pat_parsmod,
"Shear modulus (23)");
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim> void MaterialElasticOrthotropic<Dim>::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
updateInternalParameters();
AKANTU_DEBUG_OUT();
}
-/* -------------------------------------------------------------------------- */
-inline Real vector_norm(Vector<Real> & vec) {
- Real norm = 0;
- for (UInt i = 0; i < vec.size(); ++i) {
- norm += vec(i) * vec(i);
- }
- return std::sqrt(norm);
-}
-
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void MaterialElasticOrthotropic<Dim>::updateInternalParameters() {
/* 1) construction of temporary material frame stiffness tensor------------ */
// http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
Real nu21 = nu12 * E2 / E1;
Real nu31 = nu13 * E3 / E1;
Real nu32 = nu23 * E3 / E2;
// Full (i.e. dim^2 by dim^2) stiffness tensor in material frame
if (Dim == 1) {
AKANTU_DEBUG_ERROR("Dimensions 1 not implemented: makes no sense to have "
"orthotropy for 1D");
}
Real Gamma;
if (Dim == 3)
Gamma = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 -
2 * nu21 * nu32 * nu13);
if (Dim == 2)
Gamma = 1 / (1 - nu12 * nu21);
// Lamé's first parameters
this->Cprime(0, 0) = E1 * (1 - nu23 * nu32) * Gamma;
this->Cprime(1, 1) = E2 * (1 - nu13 * nu31) * Gamma;
if (Dim == 3)
this->Cprime(2, 2) = E3 * (1 - nu12 * nu21) * Gamma;
// normalised poisson's ratio's
this->Cprime(1, 0) = this->Cprime(0, 1) = E1 * (nu21 + nu31 * nu23) * Gamma;
if (Dim == 3) {
this->Cprime(2, 0) = this->Cprime(0, 2) = E1 * (nu31 + nu21 * nu32) * Gamma;
this->Cprime(2, 1) = this->Cprime(1, 2) = E2 * (nu32 + nu12 * nu31) * Gamma;
}
// Lamé's second parameters (shear moduli)
if (Dim == 3) {
this->Cprime(3, 3) = G23;
this->Cprime(4, 4) = G13;
this->Cprime(5, 5) = G12;
} else
this->Cprime(2, 2) = G12;
/* 1) rotation of C into the global frame */
this->rotateCprime();
this->C.eig(this->eigC);
}
-/* -------------------------------------------------------------------------- */
-
-template <UInt dim>
-inline void MaterialElasticOrthotropic<dim>::computePotentialEnergyOnQuad(
- const Matrix<Real> & grad_u, const Matrix<Real> & sigma, Real & epot) {
- epot = .5 * sigma.doubleDot(grad_u);
-}
-
-/* -------------------------------------------------------------------------- */
-template <UInt spatial_dimension>
-void MaterialElasticOrthotropic<spatial_dimension>::computePotentialEnergy(
- ElementType el_type, GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- Material::computePotentialEnergy(el_type, ghost_type);
-
- AKANTU_DEBUG_ASSERT(!this->finite_deformation,
- "finite deformation not possible in material orthotropic "
- "(TO BE IMPLEMENTED)");
-
- if (ghost_type != _not_ghost)
- return;
- Array<Real>::scalar_iterator epot =
- this->potential_energy(el_type, ghost_type).begin();
-
- MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
-
- computePotentialEnergyOnQuad(grad_u, sigma, *epot);
- ++epot;
-
- MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
-
- AKANTU_DEBUG_OUT();
-}
-
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
void MaterialElasticOrthotropic<spatial_dimension>::
computePotentialEnergyByElement(ElementType type, UInt index,
Vector<Real> & epot_on_quad_points) {
AKANTU_DEBUG_ASSERT(!this->finite_deformation,
"finite deformation not possible in material orthotropic "
"(TO BE IMPLEMENTED)");
Array<Real>::matrix_iterator gradu_it =
this->gradu(type).begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator gradu_end =
this->gradu(type).begin(spatial_dimension, spatial_dimension);
Array<Real>::matrix_iterator stress_it =
this->stress(type).begin(spatial_dimension, spatial_dimension);
UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type);
gradu_it += index * nb_quadrature_points;
gradu_end += (index + 1) * nb_quadrature_points;
stress_it += index * nb_quadrature_points;
Real * epot_quad = epot_on_quad_points.storage();
Matrix<Real> grad_u(spatial_dimension, spatial_dimension);
for (; gradu_it != gradu_end; ++gradu_it, ++stress_it, ++epot_quad) {
grad_u.copy(*gradu_it);
- computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad);
+ this->computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad);
}
}
/* -------------------------------------------------------------------------- */
INSTANTIATE_MATERIAL(elastic_orthotropic, MaterialElasticOrthotropic);
} // akantu
diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
index f5bad0fe2..60110342c 100644
--- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh
@@ -1,145 +1,136 @@
/**
* @file material_elastic_orthotropic.hh
*
* @author Till Junge <till.junge@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Sun Oct 19 2014
*
* @brief Orthotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "material_elastic_linear_anisotropic.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__
#define __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__
namespace akantu {
/**
* Orthotropic elastic material
*
* parameters in the material files :
* - n1 : direction of x-axis in material base, normalisation not necessary
* (default: {1, 0, 0})
* - n2 : direction of y-axis in material base, normalisation not necessary
* (default: {0, 1, 0})
* - n3 : direction of z-axis in material base, normalisation not necessary
* (default: {0, 0, 1})
* - rho : density (default: 0)
* - E1 : Young's modulus along n1 (default: 0)
* - E2 : Young's modulus along n2 (default: 0)
* - E3 : Young's modulus along n3 (default: 0)
* - nu12 : Poisson's ratio along 12 (default: 0)
* - nu13 : Poisson's ratio along 13 (default: 0)
* - nu23 : Poisson's ratio along 23 (default: 0)
* - G12 : Shear modulus along 12 (default: 0)
* - G13 : Shear modulus along 13 (default: 0)
* - G23 : Shear modulus along 23 (default: 0)
*/
template <UInt Dim>
class MaterialElasticOrthotropic
: public MaterialElasticLinearAnisotropic<Dim> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElasticOrthotropic(SolidMechanicsModel & model, const ID & id = "");
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial() override;
void updateInternalParameters() override;
- /// compute the elastic potential energy
- void computePotentialEnergy(ElementType el_type,
- GhostType ghost_type = _not_ghost) override;
-
void
computePotentialEnergyByElement(ElementType type, UInt index,
Vector<Real> & epot_on_quad_points) override;
-protected:
- inline void computePotentialEnergyOnQuad(const Matrix<Real> & grad_u,
- const Matrix<Real> & sigma,
- Real & epot);
-
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(E1, E1, Real);
AKANTU_GET_MACRO(E2, E2, Real);
AKANTU_GET_MACRO(E3, E3, Real);
AKANTU_GET_MACRO(Nu12, nu12, Real);
AKANTU_GET_MACRO(Nu13, nu13, Real);
AKANTU_GET_MACRO(Nu23, nu23, Real);
AKANTU_GET_MACRO(G12, G12, Real);
AKANTU_GET_MACRO(G13, G13, Real);
AKANTU_GET_MACRO(G23, G23, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the n1 young modulus
Real E1;
/// the n2 young modulus
Real E2;
/// the n3 young modulus
Real E3;
/// 12 Poisson coefficient
Real nu12;
/// 13 Poisson coefficient
Real nu13;
/// 23 Poisson coefficient
Real nu23;
/// 12 shear modulus
Real G12;
/// 13 shear modulus
Real G13;
/// 23 shear modulus
Real G23;
};
} // akantu
#endif /* __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__ */
diff --git a/test/test_fe_engine/_discrete_kirchhoff_triangle_18.msh b/test/test_fe_engine/_discrete_kirchhoff_triangle_18.msh
index 741047074..7422ecdb9 100644
--- a/test/test_fe_engine/_discrete_kirchhoff_triangle_18.msh
+++ b/test/test_fe_engine/_discrete_kirchhoff_triangle_18.msh
@@ -1,50 +1,271 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
-13
+100
1 0 0 0
2 1 0 0
3 1 1 0
4 0 1 0
-5 0.4999999999990217 0 0
-6 1 0.4999999999990217 0
-7 0.5000000000013878 1 0
-8 0 0.5000000000013878 0
-9 0.5000000000000262 0.5000000000000762 0
-10 0.2500000000004626 0.7500000000004626 0
-11 0.7500000000004626 0.7499999999996739 0
-12 0.7499999999996739 0.2499999999996739 0
-13 0.24999999999961 0.2500000000004752 0
+5 0.1111111111109082 0 0
+6 0.2222222222217143 0 0
+7 0.333333333332501 0 0
+8 0.4444444444432878 0 0
+9 0.5555555555543833 0 0
+10 0.6666666666657874 0 0
+11 0.7777777777771916 0 0
+12 0.8888888888885959 0 0
+13 1 0.1111111111109082 0
+14 1 0.2222222222217143 0
+15 1 0.333333333332501 0
+16 1 0.4444444444432878 0
+17 1 0.5555555555543833 0
+18 1 0.6666666666657874 0
+19 1 0.7777777777771916 0
+20 1 0.8888888888885959 0
+21 0.8888888888884262 1 0
+22 0.777777777777932 1 0
+23 0.6666666666675918 1 0
+24 0.5555555555572513 1 0
+25 0.4444444444462943 1 0
+26 0.3333333333347207 1 0
+27 0.2222222222231471 1 0
+28 0.1111111111115736 1 0
+29 0 0.8888888888884262 0
+30 0 0.777777777777932 0
+31 0 0.6666666666675918 0
+32 0 0.5555555555572513 0
+33 0 0.4444444444462943 0
+34 0 0.3333333333347207 0
+35 0 0.2222222222231471 0
+36 0 0.1111111111115736 0
+37 0.1111111111109822 0.1111111111114996 0
+38 0.1111111111110561 0.2222222222229879 0
+39 0.11111111111113 0.3333333333344741 0
+40 0.1111111111112039 0.4444444444459603 0
+41 0.1111111111112778 0.5555555555569326 0
+42 0.1111111111113518 0.6666666666673913 0
+43 0.1111111111114257 0.7777777777778498 0
+44 0.1111111111114996 0.8888888888884451 0
+45 0.2222222222218735 0.1111111111114257 0
+46 0.2222222222220326 0.2222222222228287 0
+47 0.2222222222221919 0.3333333333342274 0
+48 0.2222222222223511 0.4444444444456261 0
+49 0.2222222222225103 0.5555555555566142 0
+50 0.2222222222226695 0.6666666666671908 0
+51 0.2222222222228287 0.7777777777777674 0
+52 0.2222222222229879 0.888888888888464 0
+53 0.3333333333327477 0.1111111111113518 0
+54 0.3333333333329943 0.2222222222226695 0
+55 0.3333333333332409 0.3333333333339809 0
+56 0.3333333333334876 0.4444444444452921 0
+57 0.3333333333337343 0.5555555555562953 0
+58 0.3333333333339809 0.6666666666669903 0
+59 0.3333333333342274 0.7777777777776852 0
+60 0.3333333333344741 0.8888888888884829 0
+61 0.4444444444436218 0.1111111111112778 0
+62 0.444444444443956 0.2222222222225103 0
+63 0.44444444444429 0.3333333333337343 0
+64 0.4444444444446239 0.4444444444449581 0
+65 0.4444444444449581 0.5555555555559767 0
+66 0.4444444444452921 0.66666666666679 0
+67 0.4444444444456261 0.7777777777776029 0
+68 0.4444444444459603 0.8888888888885014 0
+69 0.5555555555547019 0.1111111111112039 0
+70 0.5555555555550206 0.2222222222223511 0
+71 0.5555555555553394 0.3333333333334876 0
+72 0.5555555555556581 0.444444444444624 0
+73 0.5555555555559767 0.5555555555556581 0
+74 0.5555555555562953 0.6666666666665895 0
+75 0.555555555556614 0.7777777777775207 0
+76 0.5555555555569326 0.8888888888885206 0
+77 0.6666666666659877 0.11111111111113 0
+78 0.6666666666661885 0.2222222222221919 0
+79 0.6666666666663887 0.3333333333332409 0
+80 0.6666666666665892 0.44444444444429 0
+81 0.66666666666679 0.5555555555553392 0
+82 0.6666666666669905 0.6666666666663887 0
+83 0.6666666666671908 0.7777777777774383 0
+84 0.6666666666673913 0.8888888888885393 0
+85 0.7777777777772739 0.1111111111110561 0
+86 0.7777777777773561 0.2222222222220327 0
+87 0.7777777777774385 0.3333333333329943 0
+88 0.7777777777775208 0.4444444444439559 0
+89 0.7777777777776028 0.5555555555550208 0
+90 0.7777777777776852 0.6666666666661885 0
+91 0.7777777777777674 0.7777777777773561 0
+92 0.7777777777778497 0.8888888888885581 0
+93 0.888888888888577 0.1111111111109821 0
+94 0.8888888888885583 0.2222222222218735 0
+95 0.8888888888885393 0.3333333333327477 0
+96 0.8888888888885205 0.4444444444436219 0
+97 0.8888888888885017 0.5555555555547019 0
+98 0.8888888888884829 0.666666666665988 0
+99 0.888888888888464 0.777777777777274 0
+100 0.8888888888884451 0.888888888888577 0
$EndNodes
$Elements
-28
-1 15 2 0 101 1
-2 15 2 0 102 2
-3 15 2 0 103 3
-4 15 2 0 104 4
-5 1 2 0 101 1 5
-6 1 2 0 101 5 2
-7 1 2 0 102 2 6
-8 1 2 0 102 6 3
-9 1 2 0 103 3 7
-10 1 2 0 103 7 4
-11 1 2 0 104 4 8
-12 1 2 0 104 8 1
-13 2 2 0 101 4 10 7
-14 2 2 0 101 1 13 8
-15 2 2 0 101 3 11 6
-16 2 2 0 101 2 12 5
-17 2 2 0 101 8 9 10
-18 2 2 0 101 7 10 9
-19 2 2 0 101 7 9 11
-20 2 2 0 101 8 13 9
-21 2 2 0 101 6 11 9
-22 2 2 0 101 5 9 13
-23 2 2 0 101 5 12 9
-24 2 2 0 101 6 9 12
-25 2 2 0 101 1 5 13
-26 2 2 0 101 4 8 10
-27 2 2 0 101 2 6 12
-28 2 2 0 101 3 7 11
+162
+1 2 2 1 101 1 5 36
+2 2 2 1 101 36 5 37
+3 2 2 1 101 36 37 35
+4 2 2 1 101 35 37 38
+5 2 2 1 101 35 38 34
+6 2 2 1 101 34 38 39
+7 2 2 1 101 34 39 33
+8 2 2 1 101 33 39 40
+9 2 2 1 101 33 40 32
+10 2 2 1 101 32 40 41
+11 2 2 1 101 32 41 31
+12 2 2 1 101 31 41 42
+13 2 2 1 101 31 42 30
+14 2 2 1 101 30 42 43
+15 2 2 1 101 30 43 29
+16 2 2 1 101 29 43 44
+17 2 2 1 101 29 44 4
+18 2 2 1 101 4 44 28
+19 2 2 1 101 5 6 37
+20 2 2 1 101 37 6 45
+21 2 2 1 101 37 45 38
+22 2 2 1 101 38 45 46
+23 2 2 1 101 38 46 39
+24 2 2 1 101 39 46 47
+25 2 2 1 101 39 47 40
+26 2 2 1 101 40 47 48
+27 2 2 1 101 40 48 41
+28 2 2 1 101 41 48 49
+29 2 2 1 101 41 49 42
+30 2 2 1 101 42 49 50
+31 2 2 1 101 42 50 43
+32 2 2 1 101 43 50 51
+33 2 2 1 101 43 51 44
+34 2 2 1 101 44 51 52
+35 2 2 1 101 44 52 28
+36 2 2 1 101 28 52 27
+37 2 2 1 101 6 7 45
+38 2 2 1 101 45 7 53
+39 2 2 1 101 45 53 46
+40 2 2 1 101 46 53 54
+41 2 2 1 101 46 54 47
+42 2 2 1 101 47 54 55
+43 2 2 1 101 47 55 48
+44 2 2 1 101 48 55 56
+45 2 2 1 101 48 56 49
+46 2 2 1 101 49 56 57
+47 2 2 1 101 49 57 50
+48 2 2 1 101 50 57 58
+49 2 2 1 101 50 58 51
+50 2 2 1 101 51 58 59
+51 2 2 1 101 51 59 52
+52 2 2 1 101 52 59 60
+53 2 2 1 101 52 60 27
+54 2 2 1 101 27 60 26
+55 2 2 1 101 7 8 53
+56 2 2 1 101 53 8 61
+57 2 2 1 101 53 61 54
+58 2 2 1 101 54 61 62
+59 2 2 1 101 54 62 55
+60 2 2 1 101 55 62 63
+61 2 2 1 101 55 63 56
+62 2 2 1 101 56 63 64
+63 2 2 1 101 56 64 57
+64 2 2 1 101 57 64 65
+65 2 2 1 101 57 65 58
+66 2 2 1 101 58 65 66
+67 2 2 1 101 58 66 59
+68 2 2 1 101 59 66 67
+69 2 2 1 101 59 67 60
+70 2 2 1 101 60 67 68
+71 2 2 1 101 60 68 26
+72 2 2 1 101 26 68 25
+73 2 2 1 101 8 9 61
+74 2 2 1 101 61 9 69
+75 2 2 1 101 61 69 62
+76 2 2 1 101 62 69 70
+77 2 2 1 101 62 70 63
+78 2 2 1 101 63 70 71
+79 2 2 1 101 63 71 64
+80 2 2 1 101 64 71 72
+81 2 2 1 101 64 72 65
+82 2 2 1 101 65 72 73
+83 2 2 1 101 65 73 66
+84 2 2 1 101 66 73 74
+85 2 2 1 101 66 74 67
+86 2 2 1 101 67 74 75
+87 2 2 1 101 67 75 68
+88 2 2 1 101 68 75 76
+89 2 2 1 101 68 76 25
+90 2 2 1 101 25 76 24
+91 2 2 1 101 9 10 69
+92 2 2 1 101 69 10 77
+93 2 2 1 101 69 77 70
+94 2 2 1 101 70 77 78
+95 2 2 1 101 70 78 71
+96 2 2 1 101 71 78 79
+97 2 2 1 101 71 79 72
+98 2 2 1 101 72 79 80
+99 2 2 1 101 72 80 73
+100 2 2 1 101 73 80 81
+101 2 2 1 101 73 81 74
+102 2 2 1 101 74 81 82
+103 2 2 1 101 74 82 75
+104 2 2 1 101 75 82 83
+105 2 2 1 101 75 83 76
+106 2 2 1 101 76 83 84
+107 2 2 1 101 76 84 24
+108 2 2 1 101 24 84 23
+109 2 2 1 101 10 11 77
+110 2 2 1 101 77 11 85
+111 2 2 1 101 77 85 78
+112 2 2 1 101 78 85 86
+113 2 2 1 101 78 86 79
+114 2 2 1 101 79 86 87
+115 2 2 1 101 79 87 80
+116 2 2 1 101 80 87 88
+117 2 2 1 101 80 88 81
+118 2 2 1 101 81 88 89
+119 2 2 1 101 81 89 82
+120 2 2 1 101 82 89 90
+121 2 2 1 101 82 90 83
+122 2 2 1 101 83 90 91
+123 2 2 1 101 83 91 84
+124 2 2 1 101 84 91 92
+125 2 2 1 101 84 92 23
+126 2 2 1 101 23 92 22
+127 2 2 1 101 11 12 85
+128 2 2 1 101 85 12 93
+129 2 2 1 101 85 93 86
+130 2 2 1 101 86 93 94
+131 2 2 1 101 86 94 87
+132 2 2 1 101 87 94 95
+133 2 2 1 101 87 95 88
+134 2 2 1 101 88 95 96
+135 2 2 1 101 88 96 89
+136 2 2 1 101 89 96 97
+137 2 2 1 101 89 97 90
+138 2 2 1 101 90 97 98
+139 2 2 1 101 90 98 91
+140 2 2 1 101 91 98 99
+141 2 2 1 101 91 99 92
+142 2 2 1 101 92 99 100
+143 2 2 1 101 92 100 22
+144 2 2 1 101 22 100 21
+145 2 2 1 101 12 2 93
+146 2 2 1 101 93 2 13
+147 2 2 1 101 93 13 94
+148 2 2 1 101 94 13 14
+149 2 2 1 101 94 14 95
+150 2 2 1 101 95 14 15
+151 2 2 1 101 95 15 96
+152 2 2 1 101 96 15 16
+153 2 2 1 101 96 16 97
+154 2 2 1 101 97 16 17
+155 2 2 1 101 97 17 98
+156 2 2 1 101 98 17 18
+157 2 2 1 101 98 18 99
+158 2 2 1 101 99 18 19
+159 2 2 1 101 99 19 100
+160 2 2 1 101 100 19 20
+161 2 2 1 101 100 20 21
+162 2 2 1 101 21 20 3
$EndElements
diff --git a/test/test_fe_engine/test_fe_engine_fixture.hh b/test/test_fe_engine/test_fe_engine_fixture.hh
index 822bf968c..8e2ba7b74 100644
--- a/test/test_fe_engine/test_fe_engine_fixture.hh
+++ b/test/test_fe_engine/test_fe_engine_fixture.hh
@@ -1,84 +1,83 @@
/* -------------------------------------------------------------------------- */
#include "element_class.hh"
#include "fe_engine.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "test_gtest_utils.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__
#define __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__
using namespace akantu;
/// Generic class for FEEngine tests
template <typename type_, template <ElementKind> class shape_t,
ElementKind kind = _ek_regular>
class TestFEMBaseFixture : public ::testing::Test {
public:
static constexpr const ElementType type = type_::value;
static constexpr const size_t dim = ElementClass<type>::getSpatialDimension();
using FEM = FEEngineTemplate<IntegratorGauss, shape_t, kind>;
/// Setup reads mesh corresponding to element type and initializes an FEEngine
void SetUp() override {
const auto dim = this->dim;
const auto type = this->type;
mesh = std::make_unique<Mesh>(dim);
std::stringstream meshfilename;
meshfilename << type << ".msh";
this->readMesh(meshfilename.str());
lower = mesh->getLowerBounds();
upper = mesh->getUpperBounds();
nb_element = this->mesh->getNbElement(type);
fem = std::make_unique<FEM>(*mesh, dim, "my_fem");
- fem->initShapeFunctions();
-
- nb_quadrature_points_total = fem->getNbIntegrationPoints(type) * nb_element;
+ nb_quadrature_points_total =
+ GaussIntegrationElement<type>::getNbQuadraturePoints() * nb_element;
SCOPED_TRACE(aka::to_string(type));
}
void TearDown() override {
fem.reset(nullptr);
mesh.reset(nullptr);
}
/// Should be reimplemented if further treatment of the mesh is needed
virtual void readMesh(std::string file_name) { mesh->read(file_name); }
protected:
std::unique_ptr<FEM> fem;
std::unique_ptr<Mesh> mesh;
UInt nb_element;
UInt nb_quadrature_points_total;
Vector<Real> lower;
Vector<Real> upper;
};
template <typename type_, template <ElementKind> class shape_t,
ElementKind kind>
constexpr const ElementType TestFEMBaseFixture<type_, shape_t, kind>::type;
template <typename type_, template <ElementKind> class shape_t,
ElementKind kind>
constexpr const size_t TestFEMBaseFixture<type_, shape_t, kind>::dim;
/* -------------------------------------------------------------------------- */
/// Base class for test with Lagrange FEEngine and regular elements
template <typename type_>
using TestFEMFixture = TestFEMBaseFixture<type_, ShapeLagrange>;
/* -------------------------------------------------------------------------- */
using types = gtest_list_t<TestElementTypes>;
TYPED_TEST_CASE(TestFEMFixture, types);
#endif /* __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__ */
diff --git a/test/test_fe_engine/test_fe_engine_gauss_integration.cc b/test/test_fe_engine/test_fe_engine_gauss_integration.cc
index 9e44bc23e..48fee7c9f 100644
--- a/test/test_fe_engine/test_fe_engine_gauss_integration.cc
+++ b/test/test_fe_engine/test_fe_engine_gauss_integration.cc
@@ -1,153 +1,155 @@
/**
* @file test_fe_engine_precomputation.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Mon Jul 13 2015
*
* @brief test integration on elements, this test consider that mesh is a cube
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "test_fe_engine_fixture.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
/* -------------------------------------------------------------------------- */
template <size_t t> using degree_t = std::integral_constant<size_t, t>;
/* -------------------------------------------------------------------------- */
using TestDegreeTypes = std::tuple<degree_t<0>, degree_t<1>, degree_t<2>, degree_t<3>, degree_t<4>, degree_t<5>>;
std::array<Polynomial<5>, 3> global_polys{
{{0.40062394, 0.13703225, 0.51731446, 0.87830084, 0.5410543, 0.71842292},
{0.41861835, 0.11080576, 0.49874043, 0.49077504, 0.85073835, 0.66259755},
{0.92620845, 0.7503478, 0.62962232, 0.31662719, 0.64069644, 0.30878135}}};
template <typename T>
class TestGaussIntegrationFixture
: public TestFEMFixture<std::tuple_element_t<0, T>> {
protected:
using parent = TestFEMFixture<std::tuple_element_t<0, T>>;
static constexpr size_t degree{std::tuple_element_t<1, T>::value};
public:
TestGaussIntegrationFixture() : integration_points_pos(0, parent::dim) {}
void SetUp() override {
parent::SetUp();
+ this->fem->initShapeFunctions();
+
auto integration_points =
this->fem->getIntegrator().template getIntegrationPoints <
parent::type,
degree == 0 ? 1 : degree > ();
nb_integration_points = integration_points.cols();
auto shapes_size = ElementClass<parent::type>::getShapeSize();
Array<Real> shapes(0, shapes_size);
this->fem->getShapeFunctions()
.template computeShapesOnIntegrationPoints<parent::type>(
this->mesh->getNodes(), integration_points, shapes, _not_ghost);
auto vect_size = this->nb_integration_points * this->nb_element;
integration_points_pos.resize(vect_size);
this->fem->getShapeFunctions()
.template interpolateOnIntegrationPoints<parent::type>(
this->mesh->getNodes(), integration_points_pos, this->dim, shapes);
for (size_t d = 0; d < this->dim; ++d) {
polys[d] = global_polys[d].extract(degree);
}
}
void testIntegrate() {
std::stringstream sstr;
sstr << this->type << ":" << this->degree;
SCOPED_TRACE(sstr.str().c_str());
auto vect_size = this->nb_integration_points * this->nb_element;
Array<Real> polynomial(vect_size);
size_t dim = parent::dim;
for (size_t d = 0; d < dim; ++d) {
auto poly = this->polys[d];
for (auto && pair :
zip(polynomial, make_view(this->integration_points_pos, dim))) {
auto && p = std::get<0>(pair);
auto & x = std::get<1>(pair);
p = poly(x(d));
}
auto res =
this->fem->getIntegrator()
.template integrate<parent::type, (degree == 0 ? 1 : degree)>(polynomial);
auto expect = poly.integrate(this->lower(d), this->upper(d));
for (size_t o = 0; o < dim; ++o) {
if (o == d)
continue;
expect *= this->upper(d) - this->lower(d);
}
EXPECT_NEAR(expect, res, 5e-14);
}
}
protected:
UInt nb_integration_points;
std::array<Array<Real>, parent::dim> polynomial;
Array<Real> integration_points_pos;
std::array<Polynomial<5>, 3> polys;
};
template <typename T>
constexpr size_t TestGaussIntegrationFixture<T>::degree;
/* -------------------------------------------------------------------------- */
/* Tests */
/* -------------------------------------------------------------------------- */
TYPED_TEST_CASE_P(TestGaussIntegrationFixture);
TYPED_TEST_P(TestGaussIntegrationFixture, ArbitraryOrder) {
this->testIntegrate();
}
REGISTER_TYPED_TEST_CASE_P(TestGaussIntegrationFixture, ArbitraryOrder);
using TestTypes =
gtest_list_t<tuple_split_t<50, cross_product_t<TestElementTypes, TestDegreeTypes>>>;
INSTANTIATE_TYPED_TEST_CASE_P(Split1, TestGaussIntegrationFixture, TestTypes);
using TestTypesTail =
gtest_list_t<tuple_split_tail_t<50, cross_product_t<TestElementTypes, TestDegreeTypes>>>;
INSTANTIATE_TYPED_TEST_CASE_P(Split2, TestGaussIntegrationFixture, TestTypesTail);
}
diff --git a/test/test_fe_engine/test_fe_engine_precomputation.cc b/test/test_fe_engine/test_fe_engine_precomputation.cc
index c66054cc6..149785f5b 100644
--- a/test/test_fe_engine/test_fe_engine_precomputation.cc
+++ b/test/test_fe_engine/test_fe_engine_precomputation.cc
@@ -1,112 +1,113 @@
/**
* @file test_fe_engine_precomputation.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Mon Jul 13 2015
*
* @brief test of the fem class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "pybind11_akantu.hh"
#include "test_fe_engine_fixture.hh"
/* -------------------------------------------------------------------------- */
#include <pybind11/embed.h>
#include <pybind11/numpy.h>
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace py = pybind11;
using namespace py::literals;
template <typename type_>
class TestFEMPyFixture : public TestFEMFixture<type_> {
using parent = TestFEMFixture<type_>;
public:
void SetUp() override {
parent::SetUp();
const auto & connectivities = this->mesh->getConnectivity(this->type);
const auto & nodes = this->mesh->getNodes().begin(this->dim);
coordinates = std::make_unique<Array<Real>>(
connectivities.size(), connectivities.getNbComponent() * this->dim);
for (auto && tuple :
zip(make_view(connectivities, connectivities.getNbComponent()),
make_view(*coordinates, this->dim,
connectivities.getNbComponent()))) {
const auto & conn = std::get<0>(tuple);
const auto & X = std::get<1>(tuple);
for (auto s : arange(conn.size())) {
Vector<Real>(X(s)) = Vector<Real>(nodes[conn(s)]);
}
}
}
void TearDown() override {
parent::TearDown();
coordinates.reset(nullptr);
}
protected:
std::unique_ptr<Array<Real>> coordinates;
};
TYPED_TEST_CASE(TestFEMPyFixture, types);
TYPED_TEST(TestFEMPyFixture, Precompute) {
SCOPED_TRACE(aka::to_string(this->type));
+ this->fem->initShapeFunctions();
const auto & N = this->fem->getShapeFunctions().getShapes(this->type);
const auto & B =
this->fem->getShapeFunctions().getShapesDerivatives(this->type);
const auto & j = this->fem->getIntegrator().getJacobians(this->type);
// Array<Real> ref_N(this->nb_quadrature_points_total, N.getNbComponent());
// Array<Real> ref_B(this->nb_quadrature_points_total, B.getNbComponent());
Array<Real> ref_j(this->nb_quadrature_points_total, j.getNbComponent());
auto ref_N(N);
auto ref_B(B);
py::module py_engine = py::module::import("py_engine");
auto py_shape = py_engine.attr("Shapes")(py::str(aka::to_string(this->type)));
auto kwargs = py::dict(
"N"_a = make_proxy(ref_N), "B"_a = make_proxy(ref_B),
"j"_a = make_proxy(ref_j), "X"_a = make_proxy(*this->coordinates),
"Q"_a = make_proxy(
this->fem->getIntegrationPoints(this->type)));
auto ret = py_shape.attr("precompute")(**kwargs);
auto check = [&](auto & ref_A, auto & A, const auto & id) {
SCOPED_TRACE(aka::to_string(this->type) + " " + id);
for (auto && n : zip(make_view(ref_A, ref_A.getNbComponent()),
make_view(A, A.getNbComponent()))) {
auto diff = (std::get<0>(n) - std::get<1>(n)).template norm<L_inf>();
EXPECT_NEAR(0., diff, 1e-10);
}
};
check(ref_N, N, "N");
check(ref_B, B, "B");
check(ref_j, j, "j");
}
diff --git a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc
index 5593e14c0..fc7e3b8e5 100644
--- a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc
+++ b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc
@@ -1,123 +1,126 @@
/**
* @file test_fe_engine_precomputation_structural.cc
*
* @author Lucas Frérot
*
* @date creation: Fri Feb 26 2018
*
* @brief test of the structural precomputations
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "fe_engine.hh"
#include "integrator_gauss.hh"
#include "shape_structural.hh"
#include "test_fe_engine_structural_fixture.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
// Need a special fixture for the extra normal
class TestBernoulliB3
: public TestFEMStructuralFixture<element_type_t<_bernoulli_beam_3>> {
using parent = TestFEMStructuralFixture<element_type_t<_bernoulli_beam_3>>;
public:
/// Load the mesh and provide extra normal direction
void readMesh(std::string filename) override {
parent::readMesh(filename);
auto & normals = this->mesh->registerData<Real>("extra_normal")
.alloc(0, dim, type, _not_ghost);
Vector<Real> normal = {-36. / 65, -48. / 65, 5. / 13};
normals.push_back(normal);
}
};
/* -------------------------------------------------------------------------- */
/// Type alias
using TestBernoulliB2 =
TestFEMStructuralFixture<element_type_t<_bernoulli_beam_2>>;
using TestDKT18 =
TestFEMStructuralFixture<element_type_t<_discrete_kirchhoff_triangle_18>>;
/* -------------------------------------------------------------------------- */
/// Solution for 2D rotation matrices
Matrix<Real> globalToLocalRotation(Real theta) {
auto c = std::cos(theta);
auto s = std::sin(theta);
return {{c, s, 0}, {-s, c, 0}, {0, 0, 1}};
}
/* -------------------------------------------------------------------------- */
TEST_F(TestBernoulliB2, PrecomputeRotations) {
+ this->fem->initShapeFunctions();
using ShapeStruct = ShapeStructural<_ek_structural>;
auto & shape = dynamic_cast<const ShapeStruct &>(fem->getShapeFunctions());
auto & rot = shape.getRotations(type);
Real a = std::atan(4. / 3);
std::vector<Real> angles = {a, -a, 0};
Math::setTolerance(1e-15);
for (auto && tuple : zip(make_view(rot, ndof, ndof), angles)) {
auto rotation = std::get<0>(tuple);
auto angle = std::get<1>(tuple);
auto rotation_error = (rotation - globalToLocalRotation(angle)).norm<L_2>();
EXPECT_NEAR(rotation_error, 0., Math::getTolerance());
}
}
/* -------------------------------------------------------------------------- */
TEST_F(TestBernoulliB3, PrecomputeRotations) {
+ this->fem->initShapeFunctions();
using ShapeStruct = ShapeStructural<_ek_structural>;
auto & shape = dynamic_cast<const ShapeStruct &>(fem->getShapeFunctions());
auto & rot = shape.getRotations(type);
Matrix<Real> ref = {{3. / 13, 4. / 13, 12. / 13},
{-4. / 5, 3. / 5, 0},
{-36. / 65, -48. / 65, 5. / 13}};
Matrix<Real> solution{ndof, ndof};
solution.block(ref, 0, 0);
solution.block(ref, dim, dim);
// The default tolerance is too much, really
Math::setTolerance(1e-15);
for (auto & rotation : make_view(rot, ndof, ndof)) {
auto rotation_error = (rotation - solution).norm<L_2>();
EXPECT_NEAR(rotation_error, 0., Math::getTolerance());
}
}
/* -------------------------------------------------------------------------- */
TEST_F(TestDKT18, PrecomputeRotations) {
+ this->fem->initShapeFunctions();
using ShapeStruct = ShapeStructural<_ek_structural>;
auto & shape = dynamic_cast<const ShapeStruct &>(fem->getShapeFunctions());
auto & rot = shape.getRotations(type);
for (auto & rotation : make_view(rot, ndof, ndof)) {
std::cout << rotation << "\n";
}
std::cout.flush();
}
diff --git a/test/test_fe_engine/test_gradient.cc b/test/test_fe_engine/test_gradient.cc
index edd30e6a9..4afb40172 100644
--- a/test/test_fe_engine/test_gradient.cc
+++ b/test/test_fe_engine/test_gradient.cc
@@ -1,103 +1,105 @@
/**
* @file test_gradient.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
*
* @date creation: Fri Sep 03 2010
* @date last modification: Thu Oct 15 2015
*
* @brief test of the fem class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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/>.
*
* @section DESCRIPTION
*
* This code is computing the gradient of a linear field and check that it gives
* a constant result. It also compute the gradient the coordinates of the mesh
* and check that it gives the identity
*
*/
/* -------------------------------------------------------------------------- */
#include "test_fe_engine_fixture.hh"
/* -------------------------------------------------------------------------- */
#include <cstdlib>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
TYPED_TEST(TestFEMFixture, GradientPoly) {
+ this->fem->initShapeFunctions();
Real alpha[2][3] = {{13, 23, 31}, {11, 7, 5}};
const auto dim = this->dim;
const auto type = this->type;
const auto & position = this->fem->getMesh().getNodes();
Array<Real> const_val(this->fem->getMesh().getNbNodes(), 2, "const_val");
for(auto && pair : zip(make_view(position, dim), make_view(const_val, 2))) {
auto & pos = std::get<0>(pair);
auto & const_ = std::get<1>(pair);
const_.set(0.);
for (UInt d = 0; d < dim; ++d) {
const_(0) += alpha[0][d] * pos(d);
const_(1) += alpha[1][d] * pos(d);
}
}
/// compute the gradient
Array<Real> grad_on_quad(this->nb_quadrature_points_total, 2 * dim, "grad_on_quad");
this->fem->gradientOnIntegrationPoints(const_val, grad_on_quad, 2, type);
/// check the results
for(auto && grad : make_view(grad_on_quad, 2, dim)) {
for (UInt d = 0; d < dim; ++d) {
EXPECT_NEAR(grad(0, d), alpha[0][d], 5e-13);
EXPECT_NEAR(grad(1, d), alpha[1][d], 5e-13);
}
}
}
TYPED_TEST(TestFEMFixture, GradientPositions) {
+ this->fem->initShapeFunctions();
const auto dim = this->dim;
const auto type = this->type;
UInt nb_quadrature_points = this->fem->getNbIntegrationPoints(type) * this->nb_element;
Array<Real> grad_coord_on_quad(nb_quadrature_points, dim * dim,
"grad_coord_on_quad");
const auto & position = this->mesh->getNodes();
this->fem->gradientOnIntegrationPoints(position, grad_coord_on_quad,
dim, type);
auto I = Matrix<Real>::eye(UInt(dim));
for(auto && grad : make_view(grad_coord_on_quad, dim, dim)) {
auto diff = (I - grad).template norm<L_inf>();
EXPECT_NEAR(0., diff, 2e-14);
}
}
}
diff --git a/test/test_fe_engine/test_integrate.cc b/test/test_fe_engine/test_integrate.cc
index 40200579e..74ad310ad 100644
--- a/test/test_fe_engine/test_integrate.cc
+++ b/test/test_fe_engine/test_integrate.cc
@@ -1,75 +1,77 @@
/**
* @file test_integrate.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
*
* @date creation: Fri Sep 03 2010
* @date last modification: Thu Oct 15 2015
*
* @brief test of the fem class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "test_fe_engine_fixture.hh"
/* -------------------------------------------------------------------------- */
#include <cstdlib>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
TYPED_TEST(TestFEMFixture, IntegrateConstant) {
+ this->fem->initShapeFunctions();
+
const auto type = this->type;
const auto & position = this->fem->getMesh().getNodes();
Array<Real> const_val(position.size(), 2, "const_val");
Array<Real> val_on_quad(this->nb_quadrature_points_total, 2, "val_on_quad");
Vector<Real> value{1, 2};
for (auto && const_ : make_view(const_val, 2)) {
const_ = value;
}
// interpolate function on quadrature points
this->fem->interpolateOnIntegrationPoints(const_val, val_on_quad, 2, type);
// integrate function on elements
Array<Real> int_val_on_elem(this->nb_element, 2, "int_val_on_elem");
this->fem->integrate(val_on_quad, int_val_on_elem, 2, type);
// get global integration value
Vector<Real> sum{0., 0.};
for (auto && int_ : make_view(int_val_on_elem, 2)) {
sum += int_;
}
auto diff = (value - sum).template norm<L_inf>();
EXPECT_NEAR(0, diff, 1e-14);
}
} // namespace
diff --git a/test/test_fe_engine/test_inverse_map.cc b/test/test_fe_engine/test_inverse_map.cc
index b4a0e9518..6da60b1df 100644
--- a/test/test_fe_engine/test_inverse_map.cc
+++ b/test/test_fe_engine/test_inverse_map.cc
@@ -1,69 +1,71 @@
/**
* @file test_inverse_map.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Fri Sep 03 2010
* @date last modification: Thu Oct 15 2015
*
* @brief test of the fem class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "test_fe_engine_fixture.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
namespace {
TYPED_TEST(TestFEMFixture, InverseMap) {
+ this->fem->initShapeFunctions();
+
Matrix<Real> quad = GaussIntegrationElement<TestFixture::type>::getQuadraturePoints();
const auto & position = this->fem->getMesh().getNodes();
/// get the quadrature points coordinates
Array<Real> coord_on_quad(quad.cols() * this->nb_element, this->dim,
"coord_on_quad");
this->fem->interpolateOnIntegrationPoints(position, coord_on_quad, this->dim, this->type);
Vector<Real> natural_coords(this->dim);
auto length = (this->upper - this->lower).template norm<L_inf>();
for(auto && enum_ : enumerate(make_view(coord_on_quad, this->dim, quad.cols()))) {
auto el = std::get<0>(enum_);
const auto & quads_coords = std::get<1>(enum_);
for (auto q : arange(quad.cols())) {
Vector<Real> quad_coord = quads_coords(q);
Vector<Real> ref_quad_coord = quad(q);
this->fem->inverseMap(quad_coord, el, this->type, natural_coords);
auto dis_normalized = ref_quad_coord.distance(natural_coords) / length;
EXPECT_NEAR(0., dis_normalized, 5e-12);
}
}
}
} // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
index 98b0cb329..ed3baeb2f 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_elastic_materials.cc
@@ -1,313 +1,879 @@
/* -------------------------------------------------------------------------- */
#include "material_elastic.hh"
#include "material_elastic_orthotropic.hh"
#include "solid_mechanics_model.hh"
#include "test_material_fixtures.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <type_traits>
/* -------------------------------------------------------------------------- */
using namespace akantu;
using types =
::testing::Types<Traits<MaterialElastic, 1>, Traits<MaterialElastic, 2>,
Traits<MaterialElastic, 3>,
Traits<MaterialElasticOrthotropic, 1>,
Traits<MaterialElasticOrthotropic, 2>,
Traits<MaterialElasticOrthotropic, 3>,
Traits<MaterialElasticLinearAnisotropic, 1>,
Traits<MaterialElasticLinearAnisotropic, 2>,
Traits<MaterialElasticLinearAnisotropic, 3>>;
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<3>>::testPushWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
+ Real E = 3.;
+ setParam("E", E);
+
+ Matrix<Real> eps = {{2}};
+ Matrix<Real> sigma(1, 1);
+ Real sigma_th = 2;
+ this->computeStressOnQuad(eps, sigma, sigma_th);
+
+ auto solution = E * eps(0, 0) + sigma_th;
+ ASSERT_NEAR(sigma(0, 0), solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<1>>::testEnergyDensity() {
+ Real eps = 2, sigma = 2;
+ Real epot = 0;
+ this->computePotentialEnergyOnQuad({{eps}}, {{sigma}}, epot);
+ Real solution = 2;
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
+ Real E = 2;
+ setParam("E", E);
+ Matrix<Real> tangent(1, 1);
+ this->computeTangentModuliOnQuad(tangent);
+ ASSERT_NEAR(tangent(0, 0), E, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<1>>::testPushWaveSpeed() {
+ Real E = 3., rho = 2.;
+ setParam("E", E);
+ setParam("rho", rho);
+
+ auto wave_speed = this->getPushWaveSpeed(Element());
+ auto solution = std::sqrt(E / rho);
+ ASSERT_NEAR(wave_speed, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<1>>::testShearWaveSpeed() {
+ ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
+ Real E = 1.;
+ Real nu = .3;
+ Real sigma_th = 0.3; // thermal stress
+ setParam("E", E);
+ setParam("nu", nu);
+
+ Matrix<Real> rotation_matrix = getRandomRotation2d();
+
+ auto grad_u = this->getDeviatoricStrain(1.).block(0, 0, 2, 2);
+
+ auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
+
+ Matrix<Real> sigma_rot(2, 2);
+ this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);
+
+ auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
+
+ Matrix<Real> identity(2, 2);
+ identity.eye();
+
+ Matrix<Real> sigma_expected =
+ 0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
+
+ auto diff = sigma - sigma_expected;
+ Real stress_error = diff.norm<L_inf>();
+ ASSERT_NEAR(stress_error, 0., 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<2>>::testEnergyDensity() {
+ Matrix<Real> sigma = {{1, 2}, {2, 4}};
+ Matrix<Real> eps = {{1, 0}, {0, 1}};
+ Real epot = 0;
+ Real solution = 2.5;
+ this->computePotentialEnergyOnQuad(eps, sigma, epot);
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <>
+void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
+ Real E = 1.;
+ Real nu = .3;
+ setParam("E", E);
+ setParam("nu", nu);
+ Matrix<Real> tangent(3, 3);
+
+ /* Plane Strain */
+ // clang-format off
+ Matrix<Real> solution = {
+ {1 - nu, nu, 0},
+ {nu, 1 - nu, 0},
+ {0, 0, (1 - 2 * nu) / 2},
+ };
+ // clang-format on
+ solution *= E / ((1 + nu) * (1 - 2 * nu));
+
+ this->computeTangentModuliOnQuad(tangent);
+ Real tangent_error = (tangent - solution).norm<L_2>();
+ ASSERT_NEAR(tangent_error, 0, 1e-14);
+
+ /* Plane Stress */
+ this->plane_stress = true;
+ this->updateInternalParameters();
+ // clang-format off
+ solution = {
+ {1, nu, 0},
+ {nu, 1, 0},
+ {0, 0, (1 - nu) / 2},
+ };
+ // clang-format on
+ solution *= E / (1 - nu * nu);
+
+ this->computeTangentModuliOnQuad(tangent);
+ tangent_error = (tangent - solution).norm<L_2>();
+ ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<2>>::testPushWaveSpeed() {
Real E = 1.;
Real nu = .3;
Real rho = 2;
setParam("E", E);
setParam("nu", nu);
setParam("rho", rho);
auto wave_speed = this->getPushWaveSpeed(Element());
Real K = E / (3 * (1 - 2 * nu));
Real mu = E / (2 * (1 + nu));
Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
ASSERT_NEAR(wave_speed, sol, 1e-14);
}
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<3>>::testShearWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<2>>::testShearWaveSpeed() {
Real E = 1.;
Real nu = .3;
Real rho = 2;
setParam("E", E);
setParam("nu", nu);
setParam("rho", rho);
auto wave_speed = this->getShearWaveSpeed(Element());
Real mu = E / (2 * (1 + nu));
Real sol = std::sqrt(mu / rho);
ASSERT_NEAR(wave_speed, sol, 1e-14);
}
/* -------------------------------------------------------------------------- */
template <> void FriendMaterial<MaterialElastic<3>>::testComputeStress() {
Real E = 1.;
Real nu = .3;
Real sigma_th = 0.3; // thermal stress
setParam("E", E);
setParam("nu", nu);
Matrix<Real> rotation_matrix = getRandomRotation3d();
auto grad_u = this->getDeviatoricStrain(1.);
auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
Matrix<Real> sigma_rot(3, 3);
this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);
auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
Matrix<Real> identity(3, 3);
identity.eye();
Matrix<Real> sigma_expected =
0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
auto diff = sigma - sigma_expected;
Real stress_error = diff.norm<L_inf>();
ASSERT_NEAR(stress_error, 0., 1e-14);
}
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElastic<3>>::testEnergyDensity() {
+ Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
+ Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+ Real epot = 0;
+ Real solution = 5.5;
+ this->computePotentialEnergyOnQuad(eps, sigma, epot);
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
+
/* -------------------------------------------------------------------------- */
template <>
void FriendMaterial<MaterialElastic<3>>::testComputeTangentModuli() {
Real E = 1.;
Real nu = .3;
setParam("E", E);
setParam("nu", nu);
Matrix<Real> tangent(6, 6);
// clang-format off
Matrix<Real> solution = {
{1 - nu, nu, nu, 0, 0, 0},
{nu, 1 - nu, nu, 0, 0, 0},
{nu, nu, 1 - nu, 0, 0, 0},
{0, 0, 0, (1 - 2 * nu) / 2, 0, 0},
{0, 0, 0, 0, (1 - 2 * nu) / 2, 0},
{0, 0, 0, 0, 0, (1 - 2 * nu) / 2},
};
// clang-format on
solution *= E / ((1 + nu) * (1 - 2 * nu));
this->computeTangentModuliOnQuad(tangent);
Real tangent_error = (tangent - solution).norm<L_2>();
ASSERT_NEAR(tangent_error, 0, 1e-14);
}
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<3>>::testEnergyDensity() {
- Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
- Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
- Real epot = 0;
- Real solution = 5.5;
- this->computePotentialEnergyOnQuad(eps, sigma, epot);
- ASSERT_NEAR(epot, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testComputeStress() {
- Real E = 3.;
- setParam("E", E);
-
- Matrix<Real> eps = {{2}};
- Matrix<Real> sigma(1, 1);
- Real sigma_th = 2;
- this->computeStressOnQuad(eps, sigma, sigma_th);
-
- auto solution = E * eps(0, 0) + sigma_th;
- ASSERT_NEAR(sigma(0, 0), solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testEnergyDensity() {
- Real eps = 2, sigma = 2;
- Real epot = 0;
- this->computePotentialEnergyOnQuad({{eps}}, {{sigma}}, epot);
- Real solution = 2;
- ASSERT_NEAR(epot, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <>
-void FriendMaterial<MaterialElastic<1>>::testComputeTangentModuli() {
- Real E = 2;
- setParam("E", E);
- Matrix<Real> tangent(1, 1);
- this->computeTangentModuliOnQuad(tangent);
- ASSERT_NEAR(tangent(0, 0), E, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testPushWaveSpeed() {
- Real E = 3., rho = 2.;
- setParam("E", E);
- setParam("rho", rho);
-
- auto wave_speed = this->getPushWaveSpeed(Element());
- auto solution = std::sqrt(E / rho);
- ASSERT_NEAR(wave_speed, solution, 1e-14);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<1>>::testShearWaveSpeed() {
- ASSERT_THROW(this->getShearWaveSpeed(Element()), debug::Exception);
-}
-
-/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testPushWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<3>>::testPushWaveSpeed() {
Real E = 1.;
Real nu = .3;
Real rho = 2;
setParam("E", E);
setParam("nu", nu);
setParam("rho", rho);
auto wave_speed = this->getPushWaveSpeed(Element());
Real K = E / (3 * (1 - 2 * nu));
Real mu = E / (2 * (1 + nu));
Real sol = std::sqrt((K + 4. / 3 * mu) / rho);
ASSERT_NEAR(wave_speed, sol, 1e-14);
}
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testShearWaveSpeed() {
+template <> void FriendMaterial<MaterialElastic<3>>::testShearWaveSpeed() {
Real E = 1.;
Real nu = .3;
Real rho = 2;
setParam("E", E);
setParam("nu", nu);
setParam("rho", rho);
auto wave_speed = this->getShearWaveSpeed(Element());
Real mu = E / (2 * (1 + nu));
Real sol = std::sqrt(mu / rho);
ASSERT_NEAR(wave_speed, sol, 1e-14);
}
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testComputeStress() {
- Real E = 1.;
- Real nu = .3;
- Real sigma_th = 0.3; // thermal stress
- setParam("E", E);
- setParam("nu", nu);
-
+template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeStress() {
+ UInt Dim = 2;
+ Real E1 = 1.;
+ Real E2 = 2.;
+ Real nu12 = 0.1;
+ Real G12 = 2.;
+
+ setParamNoUpdate("E1", E1);
+ setParamNoUpdate("E2", E2);
+ setParamNoUpdate("nu12", nu12);
+ setParamNoUpdate("G12", G12);
+
+ // material frame of reference is rotate by rotation_matrix starting from
+ // canonical basis
Matrix<Real> rotation_matrix = getRandomRotation2d();
- auto grad_u = this->getDeviatoricStrain(1.).block(0, 0, 2, 2);
+ // canonical basis as expressed in the material frame of reference, as
+ // required by MaterialElasticOrthotropic class (it is simply given by the
+ // columns of the rotation_matrix; the lines give the material basis expressed
+ // in the canonical frame of reference)
+ Vector<Real> v1(Dim);
+ Vector<Real> v2(Dim);
+ for (UInt i = 0; i < Dim; ++i) {
+ v1[i] = rotation_matrix(i, 0);
+ v2[i] = rotation_matrix(i, 1);
+ }
+
+ setParamNoUpdate("n1", v1.normalize());
+ setParamNoUpdate("n2", v2.normalize());
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
- auto grad_u_rot = this->applyRotation(grad_u, rotation_matrix);
+ // gradient in material frame of reference
+ auto grad_u = ( this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.) ).block(0, 0, 2, 2);
+
+// gradient in canonical basis (we need to rotate *back* to the canonical
+// basis)
+ auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
+// stress in the canonical basis
Matrix<Real> sigma_rot(2, 2);
- this->computeStressOnQuad(grad_u_rot, sigma_rot, sigma_th);
+ this->computeStressOnQuad(grad_u_rot, sigma_rot);
- auto sigma = this->reverseRotation(sigma_rot, rotation_matrix);
+// stress in the material reference (we need to apply the rotation)
+ auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
- Matrix<Real> identity(2, 2);
- identity.eye();
+// construction of Cijkl engineering tensor in the *material* frame of reference
+// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
+ Real nu21 = nu12*E2/E1;
+ Real gamma = 1 / ( 1 - nu12*nu21 );
- Matrix<Real> sigma_expected =
- 0.5 * E / (1 + nu) * (grad_u + grad_u.transpose()) + sigma_th * identity;
+ Matrix<Real> C_expected(2*Dim,2*Dim,0);
+ C_expected(0,0) = gamma * E1 ;
+ C_expected(1,1) = gamma * E2 ;
+ C_expected(2,2) = G12;
+
+ C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21;
+
+// epsilon is computed directly in the *material* frame of reference
+ Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
+// sigma_expected is computed directly in the *material* frame of reference
+ Matrix<Real> sigma_expected(Dim, Dim);
+ for ( UInt i = 0; i < Dim; ++i ) {
+ for ( UInt j = 0; j < Dim; ++j ) {
+ sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j);
+ }
+ }
+
+ sigma_expected(0,1) = sigma_expected(1,0) = C_expected(2,2) * 2 * epsilon(0,1);
+
+ // sigmas are checked in the *material* frame of reference
auto diff = sigma - sigma_expected;
Real stress_error = diff.norm<L_inf>();
- ASSERT_NEAR(stress_error, 0., 1e-14);
+ ASSERT_NEAR(stress_error, 0., 1e-13);
}
/* -------------------------------------------------------------------------- */
-template <>
-void FriendMaterial<MaterialElastic<2>>::testComputeTangentModuli() {
- Real E = 1.;
- Real nu = .3;
- setParam("E", E);
- setParam("nu", nu);
- Matrix<Real> tangent(3, 3);
+template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testEnergyDensity() {
+ Matrix<Real> sigma = {{1, 2}, {2, 4}};
+ Matrix<Real> eps = {{1, 0}, {0, 1}};
+ Real epot = 0;
+ Real solution = 2.5;
+ this->computePotentialEnergyOnQuad(eps, sigma, epot);
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
- /* Plane Strain */
- // clang-format off
- Matrix<Real> solution = {
- {1 - nu, nu, 0},
- {nu, 1 - nu, 0},
- {0, 0, (1 - 2 * nu) / 2},
- };
- // clang-format on
- solution *= E / ((1 + nu) * (1 - 2 * nu));
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<2>>::testComputeTangentModuli() {
+
+ // Note: for this test material and canonical basis coincide
+ Vector<Real> n1 = {1, 0};
+ Vector<Real> n2 = {0, 1};
+ Real E1 = 1.;
+ Real E2 = 2.;
+ Real nu12 = 0.1;
+ Real G12 = 2.;
+
+ setParamNoUpdate("n1", n1);
+ setParamNoUpdate("n2", n2);
+ setParamNoUpdate("E1", E1);
+ setParamNoUpdate("E2", E2);
+ setParamNoUpdate("nu12", nu12);
+ setParamNoUpdate("G12", G12);
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+// construction of Cijkl engineering tensor in the *material* frame of reference
+// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
+ Real nu21 = nu12*E2/E1;
+ Real gamma = 1 / ( 1 - nu12*nu21 );
+
+ Matrix<Real> C_expected(3,3);
+ C_expected(0,0) = gamma * E1 ;
+ C_expected(1,1) = gamma * E2 ;
+ C_expected(2,2) = G12;
+
+ C_expected(1,0) = C_expected(0,1) = gamma * E1 * nu21;
+
+ Matrix<Real> tangent(3,3);
this->computeTangentModuliOnQuad(tangent);
- Real tangent_error = (tangent - solution).norm<L_2>();
+
+ Real tangent_error = (tangent - C_expected).norm<L_2>();
ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
- /* Plane Stress */
- this->plane_stress = true;
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeStress() {
+ UInt Dim = 3;
+ Real E1 = 1.;
+ Real E2 = 2.;
+ Real E3 = 3.;
+ Real nu12 = 0.1;
+ Real nu13 = 0.2;
+ Real nu23 = 0.3;
+ Real G12 = 2.;
+ Real G13 = 3.;
+ Real G23 = 1.;
+
+ setParamNoUpdate("E1", E1);
+ setParamNoUpdate("E2", E2);
+ setParamNoUpdate("E3", E3);
+ setParamNoUpdate("nu12", nu12);
+ setParamNoUpdate("nu13", nu13);
+ setParamNoUpdate("nu23", nu23);
+ setParamNoUpdate("G12", G12);
+ setParamNoUpdate("G13", G13);
+ setParamNoUpdate("G23", G23);
+
+ // material frame of reference is rotate by rotation_matrix starting from
+ // canonical basis
+ Matrix<Real> rotation_matrix = getRandomRotation3d();
+
+ // canonical basis as expressed in the material frame of reference, as
+ // required by MaterialElasticOrthotropic class (it is simply given by the
+ // columns of the rotation_matrix; the lines give the material basis expressed
+ // in the canonical frame of reference)
+ Vector<Real> v1(Dim);
+ Vector<Real> v2(Dim);
+ Vector<Real> v3(Dim);
+ for (UInt i = 0; i < Dim; ++i) {
+ v1[i] = rotation_matrix(i, 0);
+ v2[i] = rotation_matrix(i, 1);
+ v3[i] = rotation_matrix(i, 2);
+ }
+
+ setParamNoUpdate("n1", v1.normalize());
+ setParamNoUpdate("n2", v2.normalize());
+ setParamNoUpdate("n3", v3.normalize());
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
this->updateInternalParameters();
- // clang-format off
- solution = {
- {1, nu, 0},
- {nu, 1, 0},
- {0, 0, (1 - nu) / 2},
- };
- // clang-format on
- solution *= E / (1 - nu * nu);
+ // gradient in material frame of reference
+ auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);
+
+ // gradient in canonical basis (we need to rotate *back* to the canonical
+ // basis)
+ auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
+
+ // stress in the canonical basis
+ Matrix<Real> sigma_rot(3, 3);
+ this->computeStressOnQuad(grad_u_rot, sigma_rot);
+
+ // stress in the material reference (we need to apply the rotation)
+ auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
+
+// construction of Cijkl engineering tensor in the *material* frame of reference
+// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
+ Real nu21 = nu12*E2/E1;
+ Real nu31 = nu13*E3/E1;
+ Real nu32 = nu23*E3/E2;
+ Real gamma = 1 / ( 1 - nu12*nu21 - nu23*nu32 - nu31*nu13 - 2*nu21*nu32*nu13 );
+
+ Matrix<Real> C_expected(6,6);
+ C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 );
+ C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 );
+ C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 );
+
+ C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 );
+ C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 );
+ C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 );
+
+ C_expected(3,3) = G23;
+ C_expected(4,4) = G13;
+ C_expected(5,5) = G12;
+
+ // epsilon is computed directly in the *material* frame of reference
+ Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
+
+ // sigma_expected is computed directly in the *material* frame of reference
+ Matrix<Real> sigma_expected(Dim, Dim);
+ for ( UInt i = 0; i < Dim; ++i ) {
+ for ( UInt j = 0; j < Dim; ++j ) {
+ sigma_expected(i,i) += C_expected(i,j)*epsilon(j,j);
+ }
+ }
+
+ sigma_expected(0,1) = C_expected(5,5) * 2 * epsilon(0,1);
+ sigma_expected(0,2) = C_expected(4,4) * 2 * epsilon(0,2);
+ sigma_expected(1,2) = C_expected(3,3) * 2 * epsilon(1,2);
+ sigma_expected(1,0) = sigma_expected(0,1);
+ sigma_expected(2,0) = sigma_expected(0,2);
+ sigma_expected(2,1) = sigma_expected(1,2);
+
+ // sigmas are checked in the *material* frame of reference
+ auto diff = sigma - sigma_expected;
+ Real stress_error = diff.norm<L_inf>();
+ ASSERT_NEAR(stress_error, 0., 1e-13);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testEnergyDensity() {
+ Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
+ Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+ Real epot = 0;
+ Real solution = 5.5;
+ this->computePotentialEnergyOnQuad(eps, sigma, epot);
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticOrthotropic<3>>::testComputeTangentModuli() {
+
+ // Note: for this test material and canonical basis coincide
+ UInt Dim = 3;
+ Vector<Real> n1 = {1, 0, 0};
+ Vector<Real> n2 = {0, 1, 0};
+ Vector<Real> n3 = {0, 0, 1};
+ Real E1 = 1.;
+ Real E2 = 2.;
+ Real E3 = 3.;
+ Real nu12 = 0.1;
+ Real nu13 = 0.2;
+ Real nu23 = 0.3;
+ Real G12 = 2.;
+ Real G13 = 3.;
+ Real G23 = 1.;
+
+ setParamNoUpdate("n1", n1);
+ setParamNoUpdate("n2", n2);
+ setParamNoUpdate("n3", n3);
+ setParamNoUpdate("E1", E1);
+ setParamNoUpdate("E2", E2);
+ setParamNoUpdate("E3", E3);
+ setParamNoUpdate("nu12", nu12);
+ setParamNoUpdate("nu13", nu13);
+ setParamNoUpdate("nu23", nu23);
+ setParamNoUpdate("G12", G12);
+ setParamNoUpdate("G13", G13);
+ setParamNoUpdate("G23", G23);
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+
+// construction of Cijkl engineering tensor in the *material* frame of reference
+// ref: http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13
+ Real nu21 = nu12*E2/E1;
+ Real nu31 = nu13*E3/E1;
+ Real nu32 = nu23*E3/E2;
+ Real gamma = 1 / ( 1 - nu12*nu21 - nu23*nu32 - nu31*nu13 - 2*nu21*nu32*nu13 );
+
+ Matrix<Real> C_expected(2*Dim,2*Dim,0);
+ C_expected(0,0) = gamma * E1 * ( 1-nu23*nu32 );
+ C_expected(1,1) = gamma * E2 * ( 1-nu13*nu31 );
+ C_expected(2,2) = gamma * E3 * ( 1-nu12*nu21 );
+
+ C_expected(1,0) = C_expected(0,1) = gamma * E1 * ( nu21 + nu31*nu23 );
+ C_expected(2,0) = C_expected(0,2) = gamma * E1 * ( nu31 + nu21*nu32 );
+ C_expected(2,1) = C_expected(1,2) = gamma * E2 * ( nu32 + nu12*nu31 );
+
+ C_expected(3,3) = G23;
+ C_expected(4,4) = G13;
+ C_expected(5,5) = G12;
+
+ Matrix<Real> tangent(6,6);
this->computeTangentModuliOnQuad(tangent);
- tangent_error = (tangent - solution).norm<L_2>();
+
+ Real tangent_error = (tangent - C_expected).norm<L_2>();
ASSERT_NEAR(tangent_error, 0, 1e-14);
}
/* -------------------------------------------------------------------------- */
-template <> void FriendMaterial<MaterialElastic<2>>::testEnergyDensity() {
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeStress() {
+ UInt Dim = 2;
+ Matrix<Real> C = {
+ {1.0, 0.3, 0.4},
+ {0.3, 2.0, 0.1},
+ {0.4, 0.1, 1.5},
+ };
+
+ setParamNoUpdate("C11", C(0,0));
+ setParamNoUpdate("C12", C(0,1));
+ setParamNoUpdate("C13", C(0,2));
+ setParamNoUpdate("C22", C(1,1));
+ setParamNoUpdate("C23", C(1,2));
+ setParamNoUpdate("C33", C(2,2));
+
+ // material frame of reference is rotate by rotation_matrix starting from
+ // canonical basis
+ Matrix<Real> rotation_matrix = getRandomRotation2d();
+
+ // canonical basis as expressed in the material frame of reference, as
+ // required by MaterialElasticLinearAnisotropic class (it is simply given by
+ // the columns of the rotation_matrix; the lines give the material basis
+ // expressed in the canonical frame of reference)
+ Vector<Real> v1(Dim);
+ Vector<Real> v2(Dim);
+ for (UInt i = 0; i < Dim; ++i) {
+ v1[i] = rotation_matrix(i, 0);
+ v2[i] = rotation_matrix(i, 1);
+ }
+
+ setParamNoUpdate("n1", v1.normalize());
+ setParamNoUpdate("n2", v2.normalize());
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+
+ // gradient in material frame of reference
+ auto grad_u = (this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.)).block(0, 0, 2, 2);
+
+ // gradient in canonical basis (we need to rotate *back* to the canonical
+ // basis)
+ auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
+
+ // stress in the canonical basis
+ Matrix<Real> sigma_rot(2, 2);
+ this->computeStressOnQuad(grad_u_rot, sigma_rot);
+
+ // stress in the material reference (we need to apply the rotation)
+ auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
+
+ // epsilon is computed directly in the *material* frame of reference
+ Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
+ Vector<Real> epsilon_voigt(3);
+ epsilon_voigt(0) = epsilon(0,0);
+ epsilon_voigt(1) = epsilon(1,1);
+ epsilon_voigt(2) = 2 * epsilon(0,1);
+
+ // sigma_expected is computed directly in the *material* frame of reference
+ Vector<Real> sigma_voigt = C * epsilon_voigt;
+ Matrix<Real> sigma_expected(Dim, Dim);
+ sigma_expected(0,0) = sigma_voigt(0);
+ sigma_expected(1,1) = sigma_voigt(1);
+ sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(2);
+
+ // sigmas are checked in the *material* frame of reference
+ auto diff = sigma - sigma_expected;
+ Real stress_error = diff.norm<L_inf>();
+ ASSERT_NEAR(stress_error, 0., 1e-13);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testEnergyDensity() {
Matrix<Real> sigma = {{1, 2}, {2, 4}};
Matrix<Real> eps = {{1, 0}, {0, 1}};
Real epot = 0;
Real solution = 2.5;
this->computePotentialEnergyOnQuad(eps, sigma, epot);
ASSERT_NEAR(epot, solution, 1e-14);
}
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<2>>::testComputeTangentModuli() {
+
+ // Note: for this test material and canonical basis coincide
+ Matrix<Real> C = {
+ {1.0, 0.3, 0.4},
+ {0.3, 2.0, 0.1},
+ {0.4, 0.1, 1.5},
+ };
+
+ setParamNoUpdate("C11", C(0,0));
+ setParamNoUpdate("C12", C(0,1));
+ setParamNoUpdate("C13", C(0,2));
+ setParamNoUpdate("C22", C(1,1));
+ setParamNoUpdate("C23", C(1,2));
+ setParamNoUpdate("C33", C(2,2));
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+
+ Matrix<Real> tangent(3,3);
+ this->computeTangentModuliOnQuad(tangent);
+
+ Real tangent_error = (tangent - C).norm<L_2>();
+ ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeStress() {
+ UInt Dim = 3;
+ Matrix<Real> C = {
+ {1.0, 0.3, 0.4, 0.3, 0.2, 0.1},
+ {0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
+ {0.4, 0.1, 1.5, 0.1, 0.4, 0.3},
+ {0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
+ {0.2, 0.3, 0.4, 0.1, 0.9, 0.1},
+ {0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
+ };
+
+ setParamNoUpdate("C11", C(0,0));
+ setParamNoUpdate("C12", C(0,1));
+ setParamNoUpdate("C13", C(0,2));
+ setParamNoUpdate("C14", C(0,3));
+ setParamNoUpdate("C15", C(0,4));
+ setParamNoUpdate("C16", C(0,5));
+ setParamNoUpdate("C22", C(1,1));
+ setParamNoUpdate("C23", C(1,2));
+ setParamNoUpdate("C24", C(1,3));
+ setParamNoUpdate("C25", C(1,4));
+ setParamNoUpdate("C26", C(1,5));
+ setParamNoUpdate("C33", C(2,2));
+ setParamNoUpdate("C34", C(2,3));
+ setParamNoUpdate("C35", C(2,4));
+ setParamNoUpdate("C36", C(2,5));
+ setParamNoUpdate("C44", C(3,3));
+ setParamNoUpdate("C45", C(3,4));
+ setParamNoUpdate("C46", C(3,5));
+ setParamNoUpdate("C55", C(4,4));
+ setParamNoUpdate("C56", C(4,5));
+ setParamNoUpdate("C66", C(5,5));
+
+ // material frame of reference is rotate by rotation_matrix starting from
+ // canonical basis
+ Matrix<Real> rotation_matrix = getRandomRotation3d();
+
+ // canonical basis as expressed in the material frame of reference, as
+ // required by MaterialElasticLinearAnisotropic class (it is simply given by
+ // the columns of the rotation_matrix; the lines give the material basis
+ // expressed in the canonical frame of reference)
+ Vector<Real> v1(Dim);
+ Vector<Real> v2(Dim);
+ Vector<Real> v3(Dim);
+ for (UInt i = 0; i < Dim; ++i) {
+ v1[i] = rotation_matrix(i, 0);
+ v2[i] = rotation_matrix(i, 1);
+ v3[i] = rotation_matrix(i, 2);
+ }
+
+ setParamNoUpdate("n1", v1.normalize());
+ setParamNoUpdate("n2", v2.normalize());
+ setParamNoUpdate("n3", v3.normalize());
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+
+ // gradient in material frame of reference
+ auto grad_u = this->getDeviatoricStrain(2.) + this->getHydrostaticStrain(1.);
+
+ // gradient in canonical basis (we need to rotate *back* to the canonical
+ // basis)
+ auto grad_u_rot = this->reverseRotation(grad_u, rotation_matrix);
+
+ // stress in the canonical basis
+ Matrix<Real> sigma_rot(3, 3);
+ this->computeStressOnQuad(grad_u_rot, sigma_rot);
+
+ // stress in the material reference (we need to apply the rotation)
+ auto sigma = this->applyRotation(sigma_rot, rotation_matrix);
+
+ // epsilon is computed directly in the *material* frame of reference
+ Matrix<Real> epsilon = 0.5 * ( grad_u + grad_u.transpose() );
+ Vector<Real> epsilon_voigt(6);
+ epsilon_voigt(0) = epsilon(0,0);
+ epsilon_voigt(1) = epsilon(1,1);
+ epsilon_voigt(2) = epsilon(2,2);
+ epsilon_voigt(3) = 2 * epsilon(1,2);
+ epsilon_voigt(4) = 2 * epsilon(0,2);
+ epsilon_voigt(5) = 2 * epsilon(0,1);
+
+ // sigma_expected is computed directly in the *material* frame of reference
+ Vector<Real> sigma_voigt = C * epsilon_voigt;
+ Matrix<Real> sigma_expected(Dim, Dim);
+ sigma_expected(0,0) = sigma_voigt(0);
+ sigma_expected(1,1) = sigma_voigt(1);
+ sigma_expected(2,2) = sigma_voigt(2);
+ sigma_expected(1,2) = sigma_expected(2,1) = sigma_voigt(3);
+ sigma_expected(0,2) = sigma_expected(2,0) = sigma_voigt(4);
+ sigma_expected(0,1) = sigma_expected(1,0) = sigma_voigt(5);
+
+ // sigmas are checked in the *material* frame of reference
+ auto diff = sigma - sigma_expected;
+ Real stress_error = diff.norm<L_inf>();
+ ASSERT_NEAR(stress_error, 0., 1e-13);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testEnergyDensity() {
+ Matrix<Real> sigma = {{1, 2, 3}, {2, 4, 5}, {3, 5, 6}};
+ Matrix<Real> eps = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+ Real epot = 0;
+ Real solution = 5.5;
+ this->computePotentialEnergyOnQuad(eps, sigma, epot);
+ ASSERT_NEAR(epot, solution, 1e-14);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void FriendMaterial<MaterialElasticLinearAnisotropic<3>>::testComputeTangentModuli() {
+
+ // Note: for this test material and canonical basis coincide
+ Matrix<Real> C = {
+ {1.0, 0.3, 0.4, 0.3, 0.2, 0.1},
+ {0.3, 2.0, 0.1, 0.2, 0.3, 0.2},
+ {0.4, 0.1, 1.5, 0.1, 0.4, 0.3},
+ {0.3, 0.2, 0.1, 2.4, 0.1, 0.4},
+ {0.2, 0.3, 0.4, 0.1, 0.9, 0.1},
+ {0.1, 0.2, 0.3, 0.4, 0.1, 1.2},
+ };
+
+ setParamNoUpdate("C11", C(0,0));
+ setParamNoUpdate("C12", C(0,1));
+ setParamNoUpdate("C13", C(0,2));
+ setParamNoUpdate("C14", C(0,3));
+ setParamNoUpdate("C15", C(0,4));
+ setParamNoUpdate("C16", C(0,5));
+ setParamNoUpdate("C22", C(1,1));
+ setParamNoUpdate("C23", C(1,2));
+ setParamNoUpdate("C24", C(1,3));
+ setParamNoUpdate("C25", C(1,4));
+ setParamNoUpdate("C26", C(1,5));
+ setParamNoUpdate("C33", C(2,2));
+ setParamNoUpdate("C34", C(2,3));
+ setParamNoUpdate("C35", C(2,4));
+ setParamNoUpdate("C36", C(2,5));
+ setParamNoUpdate("C44", C(3,3));
+ setParamNoUpdate("C45", C(3,4));
+ setParamNoUpdate("C46", C(3,5));
+ setParamNoUpdate("C55", C(4,4));
+ setParamNoUpdate("C56", C(4,5));
+ setParamNoUpdate("C66", C(5,5));
+
+ // set internal Cijkl matrix expressed in the canonical frame of reference
+ this->updateInternalParameters();
+
+ Matrix<Real> tangent(6,6);
+ this->computeTangentModuliOnQuad(tangent);
+
+ Real tangent_error = (tangent - C).norm<L_2>();
+ ASSERT_NEAR(tangent_error, 0, 1e-14);
+}
+
/* -------------------------------------------------------------------------- */
namespace {
template <typename T>
class TestElasticMaterialFixture : public ::TestMaterialFixture<T> {};
TYPED_TEST_CASE(TestElasticMaterialFixture, types);
TYPED_TEST(TestElasticMaterialFixture, ElasticComputeStress) {
this->material->testComputeStress();
}
TYPED_TEST(TestElasticMaterialFixture, ElasticEnergyDensity) {
this->material->testEnergyDensity();
}
TYPED_TEST(TestElasticMaterialFixture, ElasticComputeTangentModuli) {
this->material->testComputeTangentModuli();
}
TYPED_TEST(TestElasticMaterialFixture, ElasticComputePushWaveSpeed) {
this->material->testPushWaveSpeed();
}
TYPED_TEST(TestElasticMaterialFixture, ElasticComputeShearWaveSpeed) {
this->material->testShearWaveSpeed();
}
} // namespace
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
index eb4d93b77..c0a846674 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_fixtures.hh
@@ -1,171 +1,181 @@
/* -------------------------------------------------------------------------- */
#include "material_elastic.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
#include <random>
#include <type_traits>
/* -------------------------------------------------------------------------- */
#define TO_IMPLEMENT AKANTU_EXCEPTION("TEST TO IMPLEMENT")
using namespace akantu;
/* -------------------------------------------------------------------------- */
template <typename T> class FriendMaterial : public T {
public:
~FriendMaterial() = default;
virtual void testComputeStress() { TO_IMPLEMENT; };
virtual void testComputeTangentModuli() { TO_IMPLEMENT; };
virtual void testEnergyDensity() { TO_IMPLEMENT; };
virtual void testPushWaveSpeed() { TO_IMPLEMENT; }
virtual void testShearWaveSpeed() { TO_IMPLEMENT; }
static inline Matrix<Real> getDeviatoricStrain(Real intensity);
static inline Matrix<Real> getHydrostaticStrain(Real intensity);
static inline const Matrix<Real>
reverseRotation(Matrix<Real> mat, Matrix<Real> rotation_matrix) {
Matrix<Real> R = rotation_matrix;
Matrix<Real> m2 = mat * R;
Matrix<Real> m1 = R.transpose();
return m1 * m2;
};
static inline const Matrix<Real>
applyRotation(Matrix<Real> mat, const Matrix<Real> rotation_matrix) {
Matrix<Real> R = rotation_matrix;
Matrix<Real> m2 = mat * R.transpose();
Matrix<Real> m1 = R;
return m1 * m2;
};
static inline Matrix<Real> getRandomRotation3d();
static inline Matrix<Real> getRandomRotation2d();
using T::T;
};
/* -------------------------------------------------------------------------- */
template <typename T>
Matrix<Real> FriendMaterial<T>::getDeviatoricStrain(Real intensity) {
Matrix<Real> dev = {{0, 1, 2}, {1, 0, 3}, {2, 3, 0}};
dev *= intensity;
return dev;
}
/* -------------------------------------------------------------------------- */
template <typename T>
Matrix<Real> FriendMaterial<T>::getHydrostaticStrain(Real intensity) {
Matrix<Real> dev = {{1, 0, 0}, {0, 2, 0}, {0, 0, 3}};
dev *= intensity;
return dev;
}
/* -------------------------------------------------------------------------- */
template <typename T> Matrix<Real> FriendMaterial<T>::getRandomRotation3d() {
constexpr UInt Dim = 3;
// Seed with a real random value, if available
std::random_device rd;
std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
std::uniform_real_distribution<> dis;
Vector<Real> v1(Dim);
Vector<Real> v2(Dim);
Vector<Real> v3(Dim);
for (UInt i = 0; i < Dim; ++i) {
v1(i) = dis(gen);
v2(i) = dis(gen);
}
v1.normalize();
v2.normalize();
auto d = v1.dot(v2);
v2 -= v1 * d;
v2.normalize();
d = v1.dot(v2);
v2 -= v1 * d;
v2.normalize();
- v3.crossProduct(v2, v1);
+ v3.crossProduct(v1, v2);
+
+ Vector<Real> test_axis(3);
+ test_axis.crossProduct(v1, v2);
+ test_axis -= v3;
+ if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
+ AKANTU_DEBUG_ERROR("The axis vectors do not form a right-handed coordinate "
+ << "system. I. e., ||n1 x n2 - n3|| should be zero, but "
+ << "it is " << test_axis.norm() << ".");
+ }
Matrix<Real> mat(Dim, Dim);
for (UInt i = 0; i < Dim; ++i) {
mat(0, i) = v1[i];
mat(1, i) = v2[i];
mat(2, i) = v3[i];
}
return mat;
}
/* -------------------------------------------------------------------------- */
template <typename T> Matrix<Real> FriendMaterial<T>::getRandomRotation2d() {
constexpr UInt Dim = 2;
// Seed with a real random value, if available
std::random_device rd;
std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
std::uniform_real_distribution<> dis;
- Vector<Real> v1(Dim);
- Vector<Real> v2(Dim);
+ // v1 and v2 are such that they form a right-hand basis with prescribed v3,
+ // it's need (at least) for 2d linear elastic anisotropic and (orthotropic)
+ // materials to rotate the tangent stiffness
+
+ Vector<Real> v1(3);
+ Vector<Real> v2(3);
+ Vector<Real> v3 = {0,0,1};
for (UInt i = 0; i < Dim; ++i) {
v1(i) = dis(gen);
- v2(i) = dis(gen);
+ // v2(i) = dis(gen);
}
v1.normalize();
- v2.normalize();
-
- auto d = v1.dot(v2);
- v2 -= v1 * d;
- v2.normalize();
+ v2.crossProduct(v3,v1);
Matrix<Real> mat(Dim, Dim);
for (UInt i = 0; i < Dim; ++i) {
mat(0, i) = v1[i];
mat(1, i) = v2[i];
}
return mat;
}
/* -------------------------------------------------------------------------- */
template <typename T> class TestMaterialFixture : public ::testing::Test {
public:
using mat_class = typename T::mat_class;
void SetUp() override {
constexpr auto spatial_dimension = T::Dim;
mesh = std::make_unique<Mesh>(spatial_dimension);
model = std::make_unique<SolidMechanicsModel>(*mesh);
material = std::make_unique<friend_class>(*model);
}
void TearDown() override {
material.reset(nullptr);
model.reset(nullptr);
mesh.reset(nullptr);
}
using friend_class = FriendMaterial<mat_class>;
protected:
std::unique_ptr<Mesh> mesh;
std::unique_ptr<SolidMechanicsModel> model;
std::unique_ptr<friend_class> material;
};
/* -------------------------------------------------------------------------- */
template <template <UInt> class T, UInt _Dim> struct Traits {
using mat_class = T<_Dim>;
static constexpr UInt Dim = _Dim;
};
/* -------------------------------------------------------------------------- */
diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
index bac1f8c70..b2d1070d9 100644
--- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
+++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc
@@ -1,101 +1,99 @@
/**
* @file test_structural_mechanics_model_bernoulli_beam_2.cc
*
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @date creation: Fri Jul 15 2011
* @date last modification: Sun Oct 19 2014
*
* @brief Computation of the analytical exemple 1.1 in the TGC vol 6
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014, 2015 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 "test_structural_mechanics_model_fixture.hh"
/* -------------------------------------------------------------------------- */
#include <gtest/gtest.h>
using namespace akantu;
/* -------------------------------------------------------------------------- */
class TestStructDKT18
: public TestStructuralFixture<element_type_t<_discrete_kirchhoff_triangle_18>> {
using parent = TestStructuralFixture<element_type_t<_discrete_kirchhoff_triangle_18>>;
public:
void addMaterials() override {
mat.E = 1;
mat.t = 1;
mat.nu = 0.3;
this->model->addMaterial(mat);
}
void assignMaterials() override {
auto & materials = this->model->getElementMaterial(parent::type);
std::fill(materials.begin(), materials.end(), 0);
}
void setDirichlets() override {
- std::fill(this->model->getBlockedDOFs().begin(),
- this->model->getBlockedDOFs().end(), true);
+ this->model->getBlockedDOFs().set(true);
auto center_node = this->model->getBlockedDOFs().end(parent::ndof) - 1;
- // clang-format off
*center_node = {false, false, false, false, false, true};
- // clang-format on
this->model->getDisplacement().clear();
auto disp = ++this->model->getDisplacement().begin(parent::ndof);
// Displacement field from Batoz Vol. 2 p. 392
// with theta to beta correction from discrete Kirchhoff condition
- //clang-format off
+ // clang-format off
*disp = {40, 20, -800 , -20, 40, 0}; ++disp;
*disp = {50, 40, -1400, -40, 50, 0}; ++disp;
*disp = {10, 20, -200 , -20, 10, 0}; ++disp;
- //clang-format on
+ // clang-format on
}
void setNeumanns() override {}
protected:
StructuralMaterial mat;
};
/* -------------------------------------------------------------------------- */
// Batoz Vol 2. patch test, ISBN 2-86601-259-3
-
TEST_F(TestStructDKT18, TestDisplacements) {
model->solveStep();
Vector<Real> solution = {22.5, 22.5, -337.5, -22.5, 22.5, 0};
+ auto nb_nodes = this->model->getDisplacement().size();
Vector<Real> center_node_disp =
- model->getDisplacement().end(model->getSpatialDimension())[-1];
+ model->getDisplacement().begin(solution.size())[nb_nodes - 1];
auto error = solution - center_node_disp;
+ std::cout << center_node_disp << std::endl;
Real tol = Math::getTolerance();
EXPECT_NEAR(error.norm<L_2>(), 0., tol);
}