GMatElastoPlasticQPot 0.18.3
Loading...
Searching...
No Matches
Cartesian3d.h
Go to the documentation of this file.
1
7#ifndef GMATELASTOPLASTICQPOT_CARTESIAN3D_H
8#define GMATELASTOPLASTICQPOT_CARTESIAN3D_H
9
13// use "M_PI" from "math.h"
14#define _USE_MATH_DEFINES
19#include <algorithm>
20#include <math.h>
21
22#include <GMatElastic/Cartesian3d.h>
23#include <GMatTensor/Cartesian3d.h>
24#include <QPot.h>
25
26#include "config.h"
27#include "version.h"
28
29namespace GMatElastoPlasticQPot {
30
36namespace Cartesian3d {
37
48template <class T>
49inline auto Epsd(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
50{
51 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
52 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 3);
53 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 3);
54
55 using return_type = typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type;
56 return_type ret = GMatTensor::Cartesian3d::Norm_deviatoric(A);
57 ret *= std::sqrt(0.5);
58 return ret;
59}
60
67template <class T, class U>
68inline void epsd(const T& A, U& ret)
69{
70 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
71 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 3);
72 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 3);
73 GMATELASTOPLASTICQPOT_ASSERT(xt::has_shape(A, ret.shape()));
74
75 GMatTensor::Cartesian3d::norm_deviatoric(A, ret);
76 ret *= std::sqrt(0.5);
77}
78
89template <class T>
90inline auto Sigd(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
91{
92 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
93 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 3);
94 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 3);
95
96 using return_type = typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type;
97 return_type ret = GMatTensor::Cartesian3d::Norm_deviatoric(A);
98 ret *= std::sqrt(2.0);
99 return ret;
100}
101
108template <class T, class U>
109inline void sigd(const T& A, U& ret)
110{
111 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
112 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 3);
113 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 3);
114 GMATELASTOPLASTICQPOT_ASSERT(xt::has_shape(A, ret.shape()));
115
116 GMatTensor::Cartesian3d::norm_deviatoric(A, ret);
117 ret *= std::sqrt(2.0);
118}
119
128template <size_t N>
129class Cusp : public GMatElastic::Cartesian3d::Elastic<N> {
130protected:
131 array_type::tensor<size_t, N> m_i;
132 array_type::tensor<double, N> m_eps;
133 array_type::tensor<double, N + 1> m_epsy;
134 size_t m_nyield;
135
136 using GMatElastic::Cartesian3d::Elastic<N>::m_K;
137 using GMatElastic::Cartesian3d::Elastic<N>::m_G;
138 using GMatElastic::Cartesian3d::Elastic<N>::m_Eps;
139 using GMatElastic::Cartesian3d::Elastic<N>::m_Sig;
140 using GMatElastic::Cartesian3d::Elastic<N>::m_C;
141
142 using GMatTensor::Cartesian3d::Array<N>::m_ndim;
143 using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor2;
144 using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor4;
145 using GMatTensor::Cartesian3d::Array<N>::m_size;
146 using GMatTensor::Cartesian3d::Array<N>::m_shape;
147 using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor2;
148 using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor4;
149
150public:
151 using GMatTensor::Cartesian3d::Array<N>::rank;
152
153 Cusp() = default;
154
161 template <class T, class Y>
162 Cusp(const T& K, const T& G, const Y& epsy)
163 {
164 this->init_Cusp(K, G, epsy);
165 }
166
167protected:
174 template <class T, class Y>
175 void init_Cusp(const T& K, const T& G, const Y& epsy)
176 {
177 this->init_Elastic(K, G);
178 m_eps = xt::zeros<double>(m_shape);
179 m_i = xt::zeros<size_t>(m_shape);
180 this->set_epsy(epsy);
181 }
182
183public:
188 const array_type::tensor<double, N + 1>& epsy() const
189 {
190 return m_epsy;
191 }
192
197 template <class T>
198 void set_epsy(const T& epsy)
199 {
200
201 GMATELASTOPLASTICQPOT_ASSERT(epsy.dimension() == N + 1);
202 GMATELASTOPLASTICQPOT_ASSERT(
203 std::equal(m_shape.cbegin(), m_shape.cend(), epsy.shape().cbegin()));
204
205 m_epsy = epsy;
206 m_nyield = m_epsy.shape(N);
207
208#ifdef GMATELASTOPLASTICQPOT_ENABLE_ASSERT
209 for (size_t i = 0; i < m_size; ++i) {
210 double* y = &m_epsy.flat(i * m_nyield);
211 GMATELASTOPLASTICQPOT_ASSERT(std::is_sorted(y, y + m_nyield));
212 }
213#endif
214
215 this->refresh();
216 }
217
223 const array_type::tensor<size_t, N>& i() const
224 {
225 return m_i;
226 }
227
232 const array_type::tensor<double, N>& eps() const
233 {
234 return m_eps;
235 }
236
242 array_type::tensor<double, N> epsy_left() const
243 {
244 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
245
246 for (size_t i = 0; i < m_size; ++i) {
247 ret.flat(i) = m_epsy.flat(i * m_nyield + m_i.flat(i));
248 }
249
250 return ret;
251 }
252
258 array_type::tensor<double, N> epsy_right() const
259 {
260 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
261
262 for (size_t i = 0; i < m_size; ++i) {
263 ret.flat(i) = m_epsy.flat(i * m_nyield + m_i.flat(i) + 1);
264 }
265
266 return ret;
267 }
268
274 array_type::tensor<double, N> epsp() const
275 {
276 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
277
278 for (size_t i = 0; i < m_size; ++i) {
279 auto* y = &m_epsy.flat(i * m_nyield + m_i.flat(i));
280 ret.flat(i) = 0.5 * (*(y) + *(y + 1));
281 }
282
283 return ret;
284 }
285
286 void refresh(bool compute_tangent = true) override
287 {
288 (void)(compute_tangent);
289
290 namespace GT = GMatTensor::Cartesian3d::pointer;
291
292#pragma omp parallel for
293 for (size_t i = 0; i < m_size; ++i) {
294
295 double K = m_K.flat(i);
296 double G = m_G.flat(i);
297
298 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
299 double* Sig = &m_Sig.flat(i * m_stride_tensor2);
300
301 std::array<double, m_stride_tensor2> Epsd;
302 double epsm = GT::Hydrostatic_deviatoric(Eps, &Epsd[0]);
303 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&Epsd[0], &Epsd[0]));
304 m_eps.flat(i) = epsd;
305
306 const double* y = &m_epsy.flat(i * m_nyield);
307 size_t idx = QPot::iterator::lower_bound(y, y + m_nyield, epsd, m_i.flat(i));
308 m_i.flat(i) = idx;
309
310 Sig[0] = Sig[4] = Sig[8] = 3.0 * K * epsm;
311
312 if (epsd <= 0.0) {
313 Sig[1] = Sig[2] = Sig[3] = Sig[5] = Sig[6] = Sig[7] = 0.0;
314 continue;
315 }
316
317 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
318
319 double g = 2.0 * G * (1.0 - eps_min / epsd);
320 Sig[0] += g * Epsd[0];
321 Sig[1] = g * Epsd[1];
322 Sig[2] = g * Epsd[2];
323 Sig[3] = g * Epsd[3];
324 Sig[4] += g * Epsd[4];
325 Sig[5] = g * Epsd[5];
326 Sig[6] = g * Epsd[6];
327 Sig[7] = g * Epsd[7];
328 Sig[8] += g * Epsd[8];
329 }
330 }
331
332 array_type::tensor<double, N> energy() const override
333 {
334 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
335 namespace GT = GMatTensor::Cartesian3d::pointer;
336
337#pragma omp parallel for
338 for (size_t i = 0; i < m_size; ++i) {
339
340 double K = m_K.flat(i);
341 double G = m_G.flat(i);
342
343 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
344
345 std::array<double, m_stride_tensor2> Epsd;
346 double epsm = GT::Hydrostatic_deviatoric(Eps, &Epsd[0]);
347 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&Epsd[0], &Epsd[0]));
348
349 const double* y = &m_epsy.flat(i * m_nyield);
350 size_t idx = m_i.flat(i);
351
352 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
353 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
354
355 double U = 3.0 * K * std::pow(epsm, 2.0);
356 double V = 2.0 * G * (std::pow(epsd - eps_min, 2.0) - std::pow(deps_y, 2.0));
357
358 ret.flat(i) = U + V;
359 }
360
361 return ret;
362 }
363};
364
373template <size_t N>
374class Smooth : public Cusp<N> {
375protected:
376 using Cusp<N>::m_i;
377 using Cusp<N>::m_eps;
378 using Cusp<N>::m_epsy;
379 using Cusp<N>::m_nyield;
380
381 using GMatElastic::Cartesian3d::Elastic<N>::m_K;
382 using GMatElastic::Cartesian3d::Elastic<N>::m_G;
383 using GMatElastic::Cartesian3d::Elastic<N>::m_Eps;
384 using GMatElastic::Cartesian3d::Elastic<N>::m_Sig;
385 using GMatElastic::Cartesian3d::Elastic<N>::m_C;
386
387 using GMatTensor::Cartesian3d::Array<N>::m_ndim;
388 using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor2;
389 using GMatTensor::Cartesian3d::Array<N>::m_stride_tensor4;
390 using GMatTensor::Cartesian3d::Array<N>::m_size;
391 using GMatTensor::Cartesian3d::Array<N>::m_shape;
392 using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor2;
393 using GMatTensor::Cartesian3d::Array<N>::m_shape_tensor4;
394
395public:
396 using GMatTensor::Cartesian3d::Array<N>::rank;
397
398 Smooth() = default;
399
406 template <class T, class Y>
407 Smooth(const T& K, const T& G, const Y& epsy)
408 {
409 this->init_Cusp(K, G, epsy);
410 }
411
412 void refresh(bool compute_tangent = true) override
413 {
414 (void)(compute_tangent);
415
416 namespace GT = GMatTensor::Cartesian3d::pointer;
417
418#pragma omp parallel for
419 for (size_t i = 0; i < m_size; ++i) {
420
421 double K = m_K.flat(i);
422 double G = m_G.flat(i);
423
424 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
425 double* Sig = &m_Sig.flat(i * m_stride_tensor2);
426
427 std::array<double, m_stride_tensor2> Epsd;
428 double epsm = GT::Hydrostatic_deviatoric(Eps, &Epsd[0]);
429 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&Epsd[0], &Epsd[0]));
430 m_eps.flat(i) = epsd;
431
432 const double* y = &m_epsy.flat(i * m_nyield);
433 size_t idx = QPot::iterator::lower_bound(y, y + m_nyield, epsd, m_i.flat(i));
434 m_i.flat(i) = idx;
435
436 Sig[0] = Sig[4] = Sig[8] = 3.0 * K * epsm;
437
438 if (epsd <= 0.0) {
439 Sig[1] = Sig[2] = Sig[3] = Sig[5] = Sig[6] = Sig[7] = 0.0;
440 continue;
441 }
442
443 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
444 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
445
446 double g = (2.0 * G / epsd) * (deps_y / M_PI) * sin(M_PI / deps_y * (epsd - eps_min));
447 Sig[0] += g * Epsd[0];
448 Sig[1] = g * Epsd[1];
449 Sig[2] = g * Epsd[2];
450 Sig[3] = g * Epsd[3];
451 Sig[4] += g * Epsd[4];
452 Sig[5] = g * Epsd[5];
453 Sig[6] = g * Epsd[6];
454 Sig[7] = g * Epsd[7];
455 Sig[8] += g * Epsd[8];
456 }
457 }
458
459 array_type::tensor<double, N> energy() const override
460 {
461 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
462 namespace GT = GMatTensor::Cartesian3d::pointer;
463
464#pragma omp parallel for
465 for (size_t i = 0; i < m_size; ++i) {
466
467 double K = m_K.flat(i);
468 double G = m_G.flat(i);
469
470 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
471
472 std::array<double, m_stride_tensor2> Epsd;
473 double epsm = GT::Hydrostatic_deviatoric(Eps, &Epsd[0]);
474 double epsd = std::sqrt(0.5 * GT::A2s_ddot_B2s(&Epsd[0], &Epsd[0]));
475
476 const double* y = &m_epsy.flat(i * m_nyield);
477 size_t idx = m_i.flat(i);
478
479 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
480 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
481
482 double U = 3.0 * K * std::pow(epsm, 2.0);
483 double V = -4.0 * G * std::pow(deps_y / M_PI, 2.0) *
484 (1.0 + cos(M_PI / deps_y * (epsd - eps_min)));
485
486 ret.flat(i) = U + V;
487 }
488
489 return ret;
490 }
491};
492
493} // namespace Cartesian3d
494} // namespace GMatElastoPlasticQPot
495
496#endif
GMatElastic::Cartesian3d::Elastic
Array of material points with a linear elastic constitutive response.
Definition: Cartesian3d.h:103
GMatElastic::Cartesian3d::Elastic::G
const array_type::tensor< double, N > & G() const
Shear modulus per item.
Definition: Cartesian3d.h:187
GMatElastic::Cartesian3d::Elastic::Eps
const array_type::tensor< double, N+2 > & Eps() const
Strain tensor per item.
Definition: Cartesian3d.h:273
GMatElastic::Cartesian3d::Elastic::K
const array_type::tensor< double, N > & K() const
Bulk modulus per item.
Definition: Cartesian3d.h:178
GMatElastic::Cartesian3d::Elastic::m_C
array_type::tensor< double, N+4 > m_C
Tangent per item.
Definition: Cartesian3d.h:109
GMatElastic::Cartesian3d::Elastic::init_Elastic
void init_Elastic(const T &K, const T &G)
Construct system.
Definition: Cartesian3d.h:142
GMatElastic::Cartesian3d::Elastic::Sig
const array_type::tensor< double, N+2 > & Sig() const
Stress tensor per item.
Definition: Cartesian3d.h:292
GMatElastic::Cartesian3d::Elastic::m_G
array_type::tensor< double, N > m_G
Shear modulus per item.
Definition: Cartesian3d.h:106
GMatElastic::Cartesian3d::Elastic::m_Sig
array_type::tensor< double, N+2 > m_Sig
Stress tensor per item.
Definition: Cartesian3d.h:108
GMatElastic::Cartesian3d::Elastic::m_K
array_type::tensor< double, N > m_K
Bulk modulus per item.
Definition: Cartesian3d.h:105
GMatElastic::Cartesian3d::Elastic::m_Eps
array_type::tensor< double, N+2 > m_Eps
Strain tensor per item.
Definition: Cartesian3d.h:107
GMatElastoPlasticQPot::Cartesian3d::Cusp
Array of material points with an elasto-plastic material model.
Definition: Cartesian3d.h:129
GMatElastoPlasticQPot::Cartesian3d::Cusp::i
const array_type::tensor< size_t, N > & i() const
Index of the current yield strain per item.
Definition: Cartesian3d.h:223
GMatElastoPlasticQPot::Cartesian3d::Cusp::Cusp
Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian3d.h:162
GMatElastoPlasticQPot::Cartesian3d::Cusp::m_i
array_type::tensor< size_t, N > m_i
Index of the current yield strain per item.
Definition: Cartesian3d.h:131
GMatElastoPlasticQPot::Cartesian3d::Cusp::m_eps
array_type::tensor< double, N > m_eps
Equivalent strain.
Definition: Cartesian3d.h:132
GMatElastoPlasticQPot::Cartesian3d::Cusp::energy
array_type::tensor< double, N > energy() const override
Potential energy per item.
Definition: Cartesian3d.h:332
GMatElastoPlasticQPot::Cartesian3d::Cusp::init_Cusp
void init_Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian3d.h:175
GMatElastoPlasticQPot::Cartesian3d::Cusp::set_epsy
void set_epsy(const T &epsy)
Overwrite yield strains per item.
Definition: Cartesian3d.h:198
GMatElastoPlasticQPot::Cartesian3d::Cusp::epsy_right
array_type::tensor< double, N > epsy_right() const
Current yield strain right per item.
Definition: Cartesian3d.h:258
GMatElastoPlasticQPot::Cartesian3d::Cusp::epsy_left
array_type::tensor< double, N > epsy_left() const
Current yield strain left per item.
Definition: Cartesian3d.h:242
GMatElastoPlasticQPot::Cartesian3d::Cusp::epsy
const array_type::tensor< double, N+1 > & epsy() const
Yield strains per item.
Definition: Cartesian3d.h:188
GMatElastoPlasticQPot::Cartesian3d::Cusp::m_epsy
array_type::tensor< double, N+1 > m_epsy
Yield strain sequence.
Definition: Cartesian3d.h:133
GMatElastoPlasticQPot::Cartesian3d::Cusp::eps
const array_type::tensor< double, N > & eps() const
Equivalent strain per item.
Definition: Cartesian3d.h:232
GMatElastoPlasticQPot::Cartesian3d::Cusp::epsp
array_type::tensor< double, N > epsp() const
Plastic strain per item.
Definition: Cartesian3d.h:274
GMatElastoPlasticQPot::Cartesian3d::Cusp::refresh
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
Definition: Cartesian3d.h:286
GMatElastoPlasticQPot::Cartesian3d::Smooth
Array of material points with an elasto-plastic material model.
Definition: Cartesian3d.h:374
GMatElastoPlasticQPot::Cartesian3d::Smooth::energy
array_type::tensor< double, N > energy() const override
Potential energy per item.
Definition: Cartesian3d.h:459
GMatElastoPlasticQPot::Cartesian3d::Smooth::refresh
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
Definition: Cartesian3d.h:412
GMatElastoPlasticQPot::Cartesian3d::Smooth::Smooth
Smooth(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian3d.h:407
GMatTensor::Cartesian3d::Array::m_shape_tensor4
std::array< size_t, N+4 > m_shape_tensor4
Shape of an array of 4th-order tensors == [m_shape, 3, 3, 3, 3].
Definition: Cartesian3d.h:712
GMatTensor::Cartesian3d::Array::m_stride_tensor4
static constexpr size_t m_stride_tensor4
Storage stride for 4th-order tensors ( ).
Definition: Cartesian3d.h:700
GMatTensor::Cartesian3d::Array::m_ndim
static constexpr size_t m_ndim
Number of dimensions of tensors.
Definition: Cartesian3d.h:694
GMatTensor::Cartesian3d::Array::m_shape
std::array< size_t, N > m_shape
Shape of the array (of scalars).
Definition: Cartesian3d.h:706
GMatTensor::Cartesian3d::Array::m_shape_tensor2
std::array< size_t, N+2 > m_shape_tensor2
Shape of an array of 2nd-order tensors == [m_shape, 3, 3].
Definition: Cartesian3d.h:709
GMatTensor::Cartesian3d::Array::rank
static constexpr std::size_t rank
Rank of the array (the actual rank is increased with the tensor-rank).
Definition: Cartesian3d.h:597
GMatTensor::Cartesian3d::Array::m_size
size_t m_size
Size of the array (of scalars) == prod(m_shape).
Definition: Cartesian3d.h:703
GMatTensor::Cartesian3d::Array::m_stride_tensor2
static constexpr size_t m_stride_tensor2
Storage stride for 2nd-order tensors ( ).
Definition: Cartesian3d.h:697
GMatElastoPlasticQPot::Cartesian3d::Sigd
auto Sigd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent stress: norm of strain deviator.
Definition: Cartesian3d.h:90
GMatElastoPlasticQPot::Cartesian3d::Epsd
auto Epsd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent strain: norm of strain deviator.
Definition: Cartesian3d.h:49
GMatElastoPlasticQPot::Cartesian3d::epsd
void epsd(const T &A, U &ret)
Same as Epsd(), but writes to externally allocated output.
Definition: Cartesian3d.h:68
GMatElastoPlasticQPot::Cartesian3d::sigd
void sigd(const T &A, U &ret)
Same as Sigd(), but writes to externally allocated output.
Definition: Cartesian3d.h:109
GMatElastoPlasticQPot::array_type::tensor
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:60
GMatElastoPlasticQPot
Material model based on a sequence of parabolic potentials.
Definition: Cartesian2d.h:28
GMatTensor::Cartesian3d::pointer
API for individual tensors with pointer-only input.
Definition: Cartesian3d.h:724
GMatTensor::Cartesian3d::Norm_deviatoric
auto Norm_deviatoric(const T &A)
Norm of the tensor's deviator:
Definition: Cartesian3d.hpp:162
GMatTensor::Cartesian3d::norm_deviatoric
void norm_deviatoric(const T &A, R &ret)
Same as Norm_deviatoric() but writes to externally allocated output.
Definition: Cartesian3d.hpp:155
GMatTensor::allocate
Helper to allocate 'output' which is of the same type of some 'input', but of a different rank.
Definition: config.h:107
GMATELASTOPLASTICQPOT_ASSERT
#define GMATELASTOPLASTICQPOT_ASSERT(expr)
All assertions are implementation as:
Definition: config.h:30