7#ifndef GMATELASTOPLASTICQPOT_CARTESIAN2D_H
8#define GMATELASTOPLASTICQPOT_CARTESIAN2D_H
14#define _USE_MATH_DEFINES
33namespace Cartesian2d {
54 ret *= std::sqrt(0.5);
64template <
class T,
class U>
65inline void epsd(
const T& A, U& ret)
73 ret *= std::sqrt(0.5);
95 ret *= std::sqrt(2.0);
105template <
class T,
class U>
106inline void sigd(
const T& A, U& ret)
114 ret *= std::sqrt(2.0);
165 std::copy(
K.shape().cbegin(),
K.shape().cend(),
m_shape.begin());
176 auto C = xt::adapt(
m_C.data(), {m_ndim, m_ndim, m_ndim, m_ndim});
183 for (
size_t i = 0; i <
m_size; ++i) {
221 std::copy(arg.cbegin(), arg.cend(),
m_Eps.begin());
246 virtual void refresh(
bool compute_tangent =
true)
248 (void)(compute_tangent);
250#pragma omp parallel for
251 for (
size_t i = 0; i <
m_size; ++i) {
253 double K =
m_K.flat(i);
254 double G =
m_G.flat(i);
314#pragma omp parallel for
315 for (
size_t i = 0; i <
m_size; ++i) {
317 double K =
m_K.flat(i);
318 double G =
m_G.flat(i);
322 std::array<double, m_stride_tensor2>
Epsd;
323 double epsm = GT::Hydrostatic_deviatoric(
Eps, &
Epsd[0]);
324 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&
Epsd[0], &
Epsd[0]));
326 double U =
K * std::pow(epsm, 2.0);
327 double V =
G * std::pow(
epsd, 2.0);
374 template <
class T,
class Y>
387 template <
class T,
class Y>
421#ifdef GMATELASTOPLASTICQPOT_ENABLE_ASSERT
493 ret.flat(
i) = 0.5 * (*(y) + *(y + 1));
499 void refresh(
bool compute_tangent =
true)
override
501 (void)(compute_tangent);
505#pragma omp parallel for
508 double K =
m_K.flat(
i);
509 double G =
m_G.flat(
i);
514 std::array<double, m_stride_tensor2>
Epsd;
515 double epsm = GT::Hydrostatic_deviatoric(
Eps, &
Epsd[0]);
516 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&
Epsd[0], &
Epsd[0]));
520 size_t idx = QPot::iterator::lower_bound(y, y +
m_nyield,
epsd,
m_i.flat(
i));
530 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
532 double g =
G * (1.0 - eps_min /
epsd);
545#pragma omp parallel for
548 double K =
m_K.flat(
i);
549 double G =
m_G.flat(
i);
553 std::array<double, m_stride_tensor2>
Epsd;
554 double epsm = GT::Hydrostatic_deviatoric(
Eps, &
Epsd[0]);
555 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&
Epsd[0], &
Epsd[0]));
558 size_t idx =
m_i.flat(
i);
560 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
561 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
563 double U =
K * std::pow(epsm, 2.0);
564 double V =
G * (std::pow(
epsd - eps_min, 2.0) - std::pow(deps_y, 2.0));
611 template <
class T,
class Y>
617 void refresh(
bool compute_tangent =
true)
override
619 (void)(compute_tangent);
623#pragma omp parallel for
626 double K =
m_K.flat(
i);
627 double G =
m_G.flat(
i);
632 std::array<double, m_stride_tensor2>
Epsd;
633 double epsm = GT::Hydrostatic_deviatoric(
Eps, &
Epsd[0]);
634 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&
Epsd[0], &
Epsd[0]));
638 size_t idx = QPot::iterator::lower_bound(y, y +
m_nyield,
epsd,
m_i.flat(
i));
648 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
649 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
651 double g = (
G /
epsd) * (deps_y / M_PI) * sin(M_PI / deps_y * (
epsd - eps_min));
664#pragma omp parallel for
667 double K =
m_K.flat(
i);
668 double G =
m_G.flat(
i);
672 std::array<double, m_stride_tensor2>
Epsd;
673 double epsm = GT::Hydrostatic_deviatoric(
Eps, &
Epsd[0]);
674 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&
Epsd[0], &
Epsd[0]));
677 size_t idx =
m_i.flat(
i);
679 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
680 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
682 double U =
K * std::pow(epsm, 2.0);
683 double V = -2.0 *
G * std::pow(deps_y / M_PI, 2.0) *
684 (1.0 + cos(M_PI / deps_y * (
epsd - eps_min)));
Array of material points with an elasto-plastic material model.
void init_Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
array_type::tensor< double, N+1 > m_epsy
Yield strain sequence.
const array_type::tensor< size_t, N > & i() const
Index of the current yield strain per item.
array_type::tensor< double, N > energy() const override
Potential energy per item.
array_type::tensor< double, N > epsp() const
Plastic strain per item.
void set_epsy(const T &epsy)
Overwrite yield strains per item.
const array_type::tensor< double, N > & eps() const
Equivalent strain per item.
const array_type::tensor< double, N+1 > & epsy() const
Yield strains per item.
array_type::tensor< size_t, N > m_i
Index of the current yield strain per item.
Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
array_type::tensor< double, N > epsy_right() const
Current yield strain right per item.
array_type::tensor< double, N > epsy_left() const
Current yield strain left per item.
array_type::tensor< double, N > m_eps
Equivalent strain.
Array of material points with a linear elastic constitutive response.
void init_Elastic(const T &K, const T &G)
Constructor alias.
array_type::tensor< double, N > m_G
Shear modulus per item.
array_type::tensor< double, N+2 > m_Sig
Stress tensor per item.
array_type::tensor< double, N > m_K
Bulk modulus per item.
const array_type::tensor< double, N+2 > & Sig() const
Stress tensor per item.
const array_type::tensor< double, N > & G() const
Shear modulus per item.
const array_type::tensor< double, N+4 > & C() const
Tangent tensor per item.
const array_type::tensor< double, N+2 > & Eps() const
Strain tensor per item.
Elastic(const T &K, const T &G)
Construct system.
virtual void refresh(bool compute_tangent=true)
Recompute stress from strain.
array_type::tensor< double, N+4 > m_C
Tangent per item.
void set_Eps(const T &arg)
Set strain tensors.
virtual array_type::tensor< double, N > energy() const
Potential energy per item.
array_type::tensor< double, N+2 > m_Eps
Strain tensor per item.
const array_type::tensor< double, N > & K() const
Bulk modulus per item.
array_type::tensor< double, N+2 > & Eps()
Strain tensor per item.
Array of material points with an elasto-plastic material model.
array_type::tensor< double, N > energy() const override
Potential energy per item.
Smooth(const T &K, const T &G, const Y &epsy)
Construct system.
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
array_type::tensor< double, N+4 > I4d() const
Array of Cartesian2d::I4d()
std::array< size_t, N+4 > m_shape_tensor4
Shape of an array of 4th-order tensors == [m_shape, 2, 2, 2, 2].
static constexpr size_t m_stride_tensor2
Storage stride for 2nd-order tensors ( ).
size_t m_size
Size of the array (of scalars) == prod(m_shape).
array_type::tensor< double, N+4 > II() const
Array of Cartesian2d::II()
static constexpr size_t m_stride_tensor4
Storage stride for 4th-order tensors ( ).
std::array< size_t, N+2 > m_shape_tensor2
Shape of an array of 2nd-order tensors == [m_shape, 2, 2].
static constexpr size_t m_ndim
Number of dimensions of tensors.
void init(const std::array< size_t, N > &shape)
Constructor 'alias'.
std::array< size_t, N > m_shape
Shape of the array (of scalars).
static constexpr std::size_t rank
Rank of the array (the actual rank is increased with the tensor-rank).
auto Epsd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent strain: norm of strain deviator.
auto Sigd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent stress: norm of strain deviator.
void sigd(const T &A, U &ret)
Same as Sigd(), but writes to externally allocated output.
void epsd(const T &A, U &ret)
Same as Epsd(), but writes to externally allocated output.
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Material model based on a sequence of parabolic potentials.
API for individual tensors with pointer-only input.
T Hydrostatic(const T *A)
See Cartesian2d::Hydrostatic()
auto Norm_deviatoric(const T &A)
Norm of the tensor's deviator:
array_type::tensor< double, 4 > II()
Result of the dyadic product of two 2nd-order identity tensors (see I2()).
void norm_deviatoric(const T &A, R &ret)
Same as Norm_deviatoric() but writes to externally allocated output.
array_type::tensor< double, 4 > I4d()
Fourth order deviatoric projection.
Helper to allocate 'output' which is of the same type of some 'input', but of a different rank.
#define GMATELASTOPLASTICQPOT_ASSERT(expr)
All assertions are implementation as: