GMatElastoPlasticQPot 0.18.3
Loading...
Searching...
No Matches
Cartesian2d.h
Go to the documentation of this file.
1
7#ifndef GMATELASTOPLASTICQPOT_CARTESIAN2D_H
8#define GMATELASTOPLASTICQPOT_CARTESIAN2D_H
9
13// use "M_PI" from "math.h"
14#define _USE_MATH_DEFINES
19#include <algorithm>
20#include <math.h>
21
23#include <QPot.h>
24
25#include "config.h"
26#include "version.h"
27
29
33namespace Cartesian2d {
34
45template <class T>
46inline auto Epsd(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
47{
48 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
49 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 2);
50 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 2);
51
52 using return_type = typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type;
54 ret *= std::sqrt(0.5);
55 return ret;
56}
57
64template <class T, class U>
65inline void epsd(const T& A, U& ret)
66{
67 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
68 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 2);
69 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 2);
70 GMATELASTOPLASTICQPOT_ASSERT(xt::has_shape(A, ret.shape()));
71
73 ret *= std::sqrt(0.5);
74}
75
86template <class T>
87inline auto Sigd(const T& A) -> typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type
88{
89 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
90 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 2);
91 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 2);
92
93 using return_type = typename GMatTensor::allocate<xt::get_rank<T>::value - 2, T>::type;
95 ret *= std::sqrt(2.0);
96 return ret;
97}
98
105template <class T, class U>
106inline void sigd(const T& A, U& ret)
107{
108 GMATELASTOPLASTICQPOT_ASSERT(A.dimension() >= 2);
109 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 1) == 2);
110 GMATELASTOPLASTICQPOT_ASSERT(A.shape(A.dimension() - 2) == 2);
111 GMATELASTOPLASTICQPOT_ASSERT(xt::has_shape(A, ret.shape()));
112
114 ret *= std::sqrt(2.0);
115}
116
121template <size_t N>
123protected:
129
137
138public:
140
141 Elastic() = default;
142
148 template <class T>
149 Elastic(const T& K, const T& G)
150 {
151 this->init_Elastic(K, G);
152 }
153
154protected:
160 template <class T>
161 void init_Elastic(const T& K, const T& G)
162 {
163 GMATELASTOPLASTICQPOT_ASSERT(K.dimension() == N);
164 GMATELASTOPLASTICQPOT_ASSERT(xt::has_shape(K, G.shape()));
165 std::copy(K.shape().cbegin(), K.shape().cend(), m_shape.begin());
166 this->init(m_shape);
167
168 m_K = K;
169 m_G = G;
170 m_Eps = xt::zeros<double>(m_shape_tensor2);
171 m_Sig = xt::zeros<double>(m_shape_tensor2);
172 m_C = xt::empty<double>(m_shape_tensor4);
173
174#pragma omp parallel
175 {
176 auto C = xt::adapt(m_C.data(), {m_ndim, m_ndim, m_ndim, m_ndim});
177 double K;
178 double G;
181
182#pragma omp for
183 for (size_t i = 0; i < m_size; ++i) {
184 C.reset_buffer(&m_C.flat(i * m_stride_tensor4), m_stride_tensor4);
185 K = m_K.flat(i);
186 G = m_G.flat(i);
187 C = 0.5 * K * II + G * I4d;
188 }
189 }
190 }
191
192public:
198 {
199 return m_K;
200 }
201
207 {
208 return m_G;
209 }
210
217 template <class T>
218 void set_Eps(const T& arg)
219 {
221 std::copy(arg.cbegin(), arg.cend(), m_Eps.begin());
222 this->refresh();
223 }
224
246 virtual void refresh(bool compute_tangent = true)
247 {
248 (void)(compute_tangent);
249
250#pragma omp parallel for
251 for (size_t i = 0; i < m_size; ++i) {
252
253 double K = m_K.flat(i);
254 double G = m_G.flat(i);
255
256 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
257 double* Sig = &m_Sig.flat(i * m_stride_tensor2);
258
260
261 Sig[0] = (K - G) * epsm + G * Eps[0];
262 Sig[1] = G * Eps[1];
263 Sig[2] = G * Eps[2];
264 Sig[3] = (K - G) * epsm + G * Eps[3];
265 }
266 }
267
273 {
274 return m_Eps;
275 }
276
283 {
284 return m_Eps;
285 }
286
292 {
293 return m_Sig;
294 }
295
301 {
302 return m_C;
303 }
304
310 {
311 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
313
314#pragma omp parallel for
315 for (size_t i = 0; i < m_size; ++i) {
316
317 double K = m_K.flat(i);
318 double G = m_G.flat(i);
319
320 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
321
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]));
325
326 double U = K * std::pow(epsm, 2.0);
327 double V = G * std::pow(epsd, 2.0);
328
329 ret.flat(i) = U + V;
330 }
331
332 return ret;
333 }
334};
335
341template <size_t N>
342class Cusp : public Elastic<N> {
343protected:
347 size_t m_nyield;
348
349 using Elastic<N>::m_K;
350 using Elastic<N>::m_G;
351 using Elastic<N>::m_Eps;
352 using Elastic<N>::m_Sig;
353 using Elastic<N>::m_C;
354
362
363public:
365
366 Cusp() = default;
367
374 template <class T, class Y>
375 Cusp(const T& K, const T& G, const Y& epsy)
376 {
377 this->init_Cusp(K, G, epsy);
378 }
379
380protected:
387 template <class T, class Y>
388 void init_Cusp(const T& K, const T& G, const Y& epsy)
389 {
390 this->init_Elastic(K, G);
391 m_eps = xt::zeros<double>(m_shape);
392 m_i = xt::zeros<size_t>(m_shape);
393 this->set_epsy(epsy);
394 }
395
396public:
402 {
403 return m_epsy;
404 }
405
410 template <class T>
411 void set_epsy(const T& epsy)
412 {
413
414 GMATELASTOPLASTICQPOT_ASSERT(epsy.dimension() == N + 1);
416 std::equal(m_shape.cbegin(), m_shape.cend(), epsy.shape().cbegin()));
417
418 m_epsy = epsy;
419 m_nyield = m_epsy.shape(N);
420
421#ifdef GMATELASTOPLASTICQPOT_ENABLE_ASSERT
422 for (size_t i = 0; i < m_size; ++i) {
423 double* y = &m_epsy.flat(i * m_nyield);
424 GMATELASTOPLASTICQPOT_ASSERT(std::is_sorted(y, y + m_nyield));
425 }
426#endif
427
428 this->refresh();
429 }
430
437 {
438 return m_i;
439 }
440
446 {
447 return m_eps;
448 }
449
456 {
457 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
458
459 for (size_t i = 0; i < m_size; ++i) {
460 ret.flat(i) = m_epsy.flat(i * m_nyield + m_i.flat(i));
461 }
462
463 return ret;
464 }
465
472 {
473 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
474
475 for (size_t i = 0; i < m_size; ++i) {
476 ret.flat(i) = m_epsy.flat(i * m_nyield + m_i.flat(i) + 1);
477 }
478
479 return ret;
480 }
481
488 {
489 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
490
491 for (size_t i = 0; i < m_size; ++i) {
492 auto* y = &m_epsy.flat(i * m_nyield + m_i.flat(i));
493 ret.flat(i) = 0.5 * (*(y) + *(y + 1));
494 }
495
496 return ret;
497 }
498
499 void refresh(bool compute_tangent = true) override
500 {
501 (void)(compute_tangent);
502
504
505#pragma omp parallel for
506 for (size_t i = 0; i < m_size; ++i) {
507
508 double K = m_K.flat(i);
509 double G = m_G.flat(i);
510
511 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
512 double* Sig = &m_Sig.flat(i * m_stride_tensor2);
513
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]));
517 m_eps.flat(i) = epsd;
518
519 const double* y = &m_epsy.flat(i * m_nyield);
520 size_t idx = QPot::iterator::lower_bound(y, y + m_nyield, epsd, m_i.flat(i));
521 m_i.flat(i) = idx;
522
523 Sig[0] = Sig[3] = K * epsm;
524
525 if (epsd <= 0.0) {
526 Sig[1] = Sig[2] = 0.0;
527 continue;
528 }
529
530 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
531
532 double g = G * (1.0 - eps_min / epsd);
533 Sig[0] += g * Epsd[0];
534 Sig[1] = g * Epsd[1];
535 Sig[2] = g * Epsd[2];
536 Sig[3] += g * Epsd[3];
537 }
538 }
539
541 {
542 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
544
545#pragma omp parallel for
546 for (size_t i = 0; i < m_size; ++i) {
547
548 double K = m_K.flat(i);
549 double G = m_G.flat(i);
550
551 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
552
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]));
556
557 const double* y = &m_epsy.flat(i * m_nyield);
558 size_t idx = m_i.flat(i);
559
560 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
561 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
562
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));
565
566 ret.flat(i) = U + V;
567 }
568
569 return ret;
570 }
571};
572
578template <size_t N>
579class Smooth : public Cusp<N> {
580protected:
581 using Cusp<N>::m_i;
582 using Cusp<N>::m_eps;
583 using Cusp<N>::m_epsy;
584 using Cusp<N>::m_nyield;
585
586 using Elastic<N>::m_K;
587 using Elastic<N>::m_G;
588 using Elastic<N>::m_Eps;
589 using Elastic<N>::m_Sig;
590 using Elastic<N>::m_C;
591
599
600public:
602
603 Smooth() = default;
604
611 template <class T, class Y>
612 Smooth(const T& K, const T& G, const Y& epsy)
613 {
614 this->init_Cusp(K, G, epsy);
615 }
616
617 void refresh(bool compute_tangent = true) override
618 {
619 (void)(compute_tangent);
620
622
623#pragma omp parallel for
624 for (size_t i = 0; i < m_size; ++i) {
625
626 double K = m_K.flat(i);
627 double G = m_G.flat(i);
628
629 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
630 double* Sig = &m_Sig.flat(i * m_stride_tensor2);
631
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]));
635 m_eps.flat(i) = epsd;
636
637 const double* y = &m_epsy.flat(i * m_nyield);
638 size_t idx = QPot::iterator::lower_bound(y, y + m_nyield, epsd, m_i.flat(i));
639 m_i.flat(i) = idx;
640
641 Sig[0] = Sig[3] = K * epsm;
642
643 if (epsd <= 0.0) {
644 Sig[1] = Sig[2] = 0.0;
645 continue;
646 }
647
648 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
649 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
650
651 double g = (G / epsd) * (deps_y / M_PI) * sin(M_PI / deps_y * (epsd - eps_min));
652 Sig[0] += g * Epsd[0];
653 Sig[1] = g * Epsd[1];
654 Sig[2] = g * Epsd[2];
655 Sig[3] += g * Epsd[3];
656 }
657 }
658
660 {
661 array_type::tensor<double, N> ret = xt::empty<double>(m_shape);
663
664#pragma omp parallel for
665 for (size_t i = 0; i < m_size; ++i) {
666
667 double K = m_K.flat(i);
668 double G = m_G.flat(i);
669
670 const double* Eps = &m_Eps.flat(i * m_stride_tensor2);
671
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]));
675
676 const double* y = &m_epsy.flat(i * m_nyield);
677 size_t idx = m_i.flat(i);
678
679 double eps_min = 0.5 * (*(y + idx) + *(y + idx + 1));
680 double deps_y = 0.5 * (*(y + idx + 1) - *(y + idx));
681
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)));
685
686 ret.flat(i) = U + V;
687 }
688
689 return ret;
690 }
691};
692
693} // namespace Cartesian2d
694} // namespace GMatElastoPlasticQPot
695
696#endif
Array of material points with an elasto-plastic material model.
Definition: Cartesian2d.h:342
void init_Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian2d.h:388
array_type::tensor< double, N+1 > m_epsy
Yield strain sequence.
Definition: Cartesian2d.h:346
const array_type::tensor< size_t, N > & i() const
Index of the current yield strain per item.
Definition: Cartesian2d.h:436
array_type::tensor< double, N > energy() const override
Potential energy per item.
Definition: Cartesian2d.h:540
array_type::tensor< double, N > epsp() const
Plastic strain per item.
Definition: Cartesian2d.h:487
void set_epsy(const T &epsy)
Overwrite yield strains per item.
Definition: Cartesian2d.h:411
const array_type::tensor< double, N > & eps() const
Equivalent strain per item.
Definition: Cartesian2d.h:445
const array_type::tensor< double, N+1 > & epsy() const
Yield strains per item.
Definition: Cartesian2d.h:401
array_type::tensor< size_t, N > m_i
Index of the current yield strain per item.
Definition: Cartesian2d.h:344
Cusp(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian2d.h:375
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
Definition: Cartesian2d.h:499
array_type::tensor< double, N > epsy_right() const
Current yield strain right per item.
Definition: Cartesian2d.h:471
array_type::tensor< double, N > epsy_left() const
Current yield strain left per item.
Definition: Cartesian2d.h:455
array_type::tensor< double, N > m_eps
Equivalent strain.
Definition: Cartesian2d.h:345
Array of material points with a linear elastic constitutive response.
Definition: Cartesian2d.h:122
void init_Elastic(const T &K, const T &G)
Constructor alias.
Definition: Cartesian2d.h:161
array_type::tensor< double, N > m_G
Shear modulus per item.
Definition: Cartesian2d.h:125
array_type::tensor< double, N+2 > m_Sig
Stress tensor per item.
Definition: Cartesian2d.h:127
array_type::tensor< double, N > m_K
Bulk modulus per item.
Definition: Cartesian2d.h:124
const array_type::tensor< double, N+2 > & Sig() const
Stress tensor per item.
Definition: Cartesian2d.h:291
const array_type::tensor< double, N > & G() const
Shear modulus per item.
Definition: Cartesian2d.h:206
const array_type::tensor< double, N+4 > & C() const
Tangent tensor per item.
Definition: Cartesian2d.h:300
const array_type::tensor< double, N+2 > & Eps() const
Strain tensor per item.
Definition: Cartesian2d.h:272
Elastic(const T &K, const T &G)
Construct system.
Definition: Cartesian2d.h:149
virtual void refresh(bool compute_tangent=true)
Recompute stress from strain.
Definition: Cartesian2d.h:246
array_type::tensor< double, N+4 > m_C
Tangent per item.
Definition: Cartesian2d.h:128
void set_Eps(const T &arg)
Set strain tensors.
Definition: Cartesian2d.h:218
virtual array_type::tensor< double, N > energy() const
Potential energy per item.
Definition: Cartesian2d.h:309
array_type::tensor< double, N+2 > m_Eps
Strain tensor per item.
Definition: Cartesian2d.h:126
const array_type::tensor< double, N > & K() const
Bulk modulus per item.
Definition: Cartesian2d.h:197
array_type::tensor< double, N+2 > & Eps()
Strain tensor per item.
Definition: Cartesian2d.h:282
Array of material points with an elasto-plastic material model.
Definition: Cartesian2d.h:579
array_type::tensor< double, N > energy() const override
Potential energy per item.
Definition: Cartesian2d.h:659
Smooth(const T &K, const T &G, const Y &epsy)
Construct system.
Definition: Cartesian2d.h:612
void refresh(bool compute_tangent=true) override
Recompute stress from strain.
Definition: Cartesian2d.h:617
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].
Definition: Cartesian2d.h:600
static constexpr size_t m_stride_tensor2
Storage stride for 2nd-order tensors ( ).
Definition: Cartesian2d.h:585
size_t m_size
Size of the array (of scalars) == prod(m_shape).
Definition: Cartesian2d.h:591
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 ( ).
Definition: Cartesian2d.h:588
std::array< size_t, N+2 > m_shape_tensor2
Shape of an array of 2nd-order tensors == [m_shape, 2, 2].
Definition: Cartesian2d.h:597
static constexpr size_t m_ndim
Number of dimensions of tensors.
Definition: Cartesian2d.h:582
void init(const std::array< size_t, N > &shape)
Constructor 'alias'.
std::array< size_t, N > m_shape
Shape of the array (of scalars).
Definition: Cartesian2d.h:594
static constexpr std::size_t rank
Rank of the array (the actual rank is increased with the tensor-rank).
Definition: Cartesian2d.h:485
auto Epsd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent strain: norm of strain deviator.
Definition: Cartesian2d.h:46
auto Sigd(const T &A) -> typename GMatTensor::allocate< xt::get_rank< T >::value - 2, T >::type
Equivalent stress: norm of strain deviator.
Definition: Cartesian2d.h:87
void sigd(const T &A, U &ret)
Same as Sigd(), but writes to externally allocated output.
Definition: Cartesian2d.h:106
void epsd(const T &A, U &ret)
Same as Epsd(), but writes to externally allocated output.
Definition: Cartesian2d.h:65
xt::xtensor< T, N > tensor
Fixed (static) rank array.
Definition: config.h:60
Material model based on a sequence of parabolic potentials.
Definition: Cartesian2d.h:28
API for individual tensors with pointer-only input.
Definition: Cartesian2d.h:612
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()).
Definition: Cartesian2d.hpp:54
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.
Definition: Cartesian2d.hpp:82
Helper to allocate 'output' which is of the same type of some 'input', but of a different rank.
Definition: config.h:107
#define GMATELASTOPLASTICQPOT_ASSERT(expr)
All assertions are implementation as:
Definition: config.h:30