MiniTensor Version of the Day
Loading...
Searching...
No Matches
MiniTensor_Mechanics.h
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MiniTensor Package
4//
5// Copyright 2016 NTESS and the MiniTensor contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#if !defined(MiniTensor_Mechanics_h)
11#define MiniTensor_Mechanics_h
12
13#include "MiniTensor_Tensor.h"
14
15namespace minitensor {
16
19
25template<typename T, Index N>
27Tensor<T, N>
28vol(Tensor<T, N> const & A);
29
35template<typename T, Index N>
37Tensor<T, N>
38dev(Tensor<T, N> const & A);
39
46template<typename T, Index N>
48Vector<T, N>
49push_forward_covariant(Tensor<T, N> const & F, Vector<T, N> const & u);
50
57template<typename T, Index N>
59Vector<T, N>
60pull_back_covariant(Tensor<T, N> const & F, Vector<T, N> const & u);
61
68template<typename T, Index N>
70Vector<T, N>
71push_forward_contravariant(Tensor<T, N> const & F, Vector<T, N> const & u);
72
79template<typename T, Index N>
81Vector<T, N>
82pull_back_contravariant(Tensor<T, N> const & F, Vector<T, N> const & u);
83
90template<typename T, Index N>
92Tensor<T, N>
93push_forward_covariant(Tensor<T, N> const & F, Tensor<T, N> const & A);
94
101template<typename T, Index N>
103Tensor<T, N>
104pull_back_covariant(Tensor<T, N> const & F, Tensor<T, N> const & A);
105
112template<typename T, Index N>
114Tensor<T, N>
115push_forward_contravariant(Tensor<T, N> const & F, Tensor<T, N> const & A);
116
123template<typename T, Index N>
125Tensor<T, N>
126pull_back_contravariant(Tensor<T, N> const & F, Tensor<T, N> const & A);
127
134template<typename T, Index N>
136Vector<T, N>
137piola(Tensor<T, N> const & F, Vector<T, N> const & u);
138
145template<typename T, Index N>
147Vector<T, N>
148piola_inverse(Tensor<T, N> const & F, Vector<T, N> const & u);
149
157template<typename T, Index N>
159Tensor<T, N>
160piola(Tensor<T, N> const & F, Tensor<T, N> const & sigma);
161
169template<typename T, Index N>
171Tensor<T, N>
172piola_inverse(Tensor<T, N> const & F, Tensor<T, N> const & P);
173
177template<typename T, Index N>
179T
180smallest_eigenvalue(Tensor<T, N> const & A);
181
190template<typename T, Index N>
192bool
193check_strict_ellipticity(Tensor4<T, N> const & A);
194
201template<typename T, Index N>
202std::pair<bool, Vector<T, N>>
203check_strong_ellipticity(Tensor4<T, N> const & A);
204
205} // namespace minitensor
206
207namespace minitensor {
208
209//
210// R^N volumetric part of 2nd-order tensor
211// \return \f$ \frac{1}{N} \mathrm{tr}\:(A) I \f$
212//
213template<typename T, Index N>
215Tensor<T, N>
217{
218 Index const
219 dimension = A.get_dimension();
220
221 T const
222 theta = (1.0/dimension) * trace(A);
223
224 return theta * eye<T, N>(dimension);
225}
226
227//
228// R^N deviatoric part of 2nd-order tensor
229// \return \f$ A - vol(A) \f$
230//
231template<typename T, Index N>
233Tensor<T, N>
235{
236 return A - vol(A);
237}
238
239} // namespace minitensor
240namespace minitensor {
241
242//
243// Push forward covariant vector
244// \param \f$ F, u \f$
245// \return \f$ F^{-T} u \f$
246//
247template<typename T, Index N>
249Vector<T, N>
251{
252 Index const
253 dimension = F.get_dimension();
254
256 v(dimension);
257
258 T const
259 J = det(F);
260
261 assert(J > 0.0);
262
263 switch (dimension) {
264
265 default:
266 MT_ERROR_EXIT("Supports only 2D and 3D.");
267 break;
268
269 case 3:
270 v(0) = (
271 (-F(1,2)*F(2,1) + F(1,1)*F(2,2)) * u(0) +
272 ( F(1,2)*F(2,0) - F(1,0)*F(2,2)) * u(1) +
273 (-F(1,1)*F(2,0) + F(1,0)*F(2,1)) * u(2)) / J;
274
275 v(1) = (
276 ( F(0,2)*F(2,1) - F(0,1)*F(2,2)) * u(0) +
277 (-F(0,2)*F(2,0) + F(0,0)*F(2,2)) * u(1) +
278 ( F(0,1)*F(2,0) - F(0,0)*F(2,1)) * u(2)) / J;
279
280 v(2) = (
281 (-F(0,2)*F(1,1) + F(0,1)*F(1,2)) * u(0) +
282 ( F(0,2)*F(1,0) - F(0,0)*F(1,2)) * u(1) +
283 (-F(0,1)*F(1,0) + F(0,0)*F(1,1)) * u(2)) / J;
284
285 break;
286
287 case 2:
288 v(0) = ( F(1,1) * u(0) - F(1,0) * u(1)) / J;
289 v(1) = (-F(0,1) * u(0) + F(0,0) * u(1)) / J;
290 break;
291
292 }
293
294 return v;
295}
296
297//
298// Pull back covariant vector
299// \param \f$ F, v \f$
300// \return \f$ F^T v \f$
301//
302template<typename T, Index N>
304Vector<T, N>
306{
307 Index const
308 dimension = F.get_dimension();
309
311 v(dimension);
312
313 switch (dimension) {
314
315 default:
316 MT_ERROR_EXIT("Supports only 2D and 3D.");
317 break;
318
319 case 3:
320 v(0) = F(0,0) * u(0) + F(1,0) * u(1) + F(2,0) * u(2);
321 v(1) = F(0,1) * u(0) + F(1,1) * u(1) + F(2,1) * u(2);
322 v(2) = F(0,2) * u(0) + F(1,2) * u(1) + F(2,2) * u(2);
323
324 break;
325
326 case 2:
327 v(0) = F(0,0) * u(0) + F(1,0) * u(1);
328 v(1) = F(0,1) * u(0) + F(1,1) * u(1);
329
330 break;
331
332 }
333
334 return v;
335}
336
337//
338// Push forward contravariant vector
339// \param \f$ F, u \f$
340// \return \f$ F u \f$
341//
342template<typename T, Index N>
344Vector<T, N>
346{
347 Index const
348 dimension = F.get_dimension();
349
351 v(dimension);
352
353 switch (dimension) {
354
355 default:
356 MT_ERROR_EXIT("Supports only 2D and 3D.");
357 break;
358
359 case 3:
360 v(0) = F(0,0) * u(0) + F(0,1) * u(1) + F(0,2) * u(2);
361 v(1) = F(1,0) * u(0) + F(1,1) * u(1) + F(1,2) * u(2);
362 v(2) = F(2,0) * u(0) + F(2,1) * u(1) + F(2,2) * u(2);
363
364 break;
365
366 case 2:
367 v(0) = F(0,0) * u(0) + F(0,1) * u(1);
368 v(1) = F(1,0) * u(0) + F(1,1) * u(1);
369
370 break;
371
372 }
373
374 return v;
375}
376
377//
378// Pull back contravariant vector
379// \param \f$ F, u \f$
380// \return \f$ F^{-1} u \f$
381//
382template<typename T, Index N>
384Vector<T, N>
386{
387 Index const
388 dimension = F.get_dimension();
389
391 v(dimension);
392
393 T const
394 J = det(F);
395
396 assert(J > 0.0);
397
398 switch (dimension) {
399
400 default:
401 MT_ERROR_EXIT("Supports only 2D and 3D.");
402 break;
403
404 case 3:
405 v(0) = (
406 (-F(1,2)*F(2,1) + F(1,1)*F(2,2)) * u(0) +
407 ( F(0,2)*F(2,1) - F(0,1)*F(2,2)) * u(1) +
408 (-F(0,2)*F(1,1) + F(0,1)*F(1,2)) * u(2)) / J;
409
410 v(1) = (
411 ( F(1,2)*F(2,0) - F(1,0)*F(2,2)) * u(0) +
412 (-F(0,2)*F(2,0) + F(0,0)*F(2,2)) * u(1) +
413 ( F(0,2)*F(1,0) - F(0,0)*F(1,2)) * u(2)) / J;
414
415 v(2) = (
416 (-F(1,1)*F(2,0) + F(1,0)*F(2,1)) * u(0) +
417 ( F(0,1)*F(2,0) - F(0,0)*F(2,1)) * u(1) +
418 (-F(0,1)*F(1,0) + F(0,0)*F(1,1)) * u(2)) / J;
419
420 break;
421
422 case 2:
423 v(0) = ( F(1,1) * u(0) - F(0,1) * u(1)) / J;
424 v(1) = (-F(1,0) * u(0) + F(0,0) * u(1)) / J;
425 break;
426
427 }
428
429 return v;
430}
431
432//
433// Push forward covariant tensor
434// \param \f$ F, A \f$
435// \return \f$ F^{-T} A F^{-1} \f$
436//
437template<typename T, Index N>
439Tensor<T, N>
441{
442 Index const
443 dimension = F.get_dimension();
444
446 G(dimension);
447
448 T const
449 J = det(F);
450
451 assert(J > 0.0);
452
453 switch (dimension) {
454
455 default:
456 MT_ERROR_EXIT("Supports only 2D and 3D.");
457 break;
458
459 case 3:
460 G(0,0) = (-F(1,2)*F(2,1) + F(1,1)*F(2,2)) / J;
461 G(0,1) = ( F(0,2)*F(2,1) - F(0,1)*F(2,2)) / J;
462 G(0,2) = (-F(0,2)*F(1,1) + F(0,1)*F(1,2)) / J;
463
464 G(1,0) = ( F(1,2)*F(2,0) - F(1,0)*F(2,2)) / J;
465 G(1,1) = (-F(0,2)*F(2,0) + F(0,0)*F(2,2)) / J;
466 G(1,2) = ( F(0,2)*F(1,0) - F(0,0)*F(1,2)) / J;
467
468 G(2,0) = (-F(1,1)*F(2,0) + F(1,0)*F(2,1)) / J;
469 G(2,1) = ( F(0,1)*F(2,0) - F(0,0)*F(2,1)) / J;
470 G(2,2) = (-F(0,1)*F(1,0) + F(0,0)*F(1,1)) / J;
471 break;
472
473 case 2:
474 G(0,0) = F(1,1) / J;
475 G(0,1) = -F(0,1) / J;
476
477 G(1,0) = -F(1,0) / J;
478 G(1,1) = F(0,0) / J;
479 break;
480
481 }
482
483 return t_dot(G, dot(A, G));
484}
485
486//
487// Pull back covariant tensor
488// \param \f$ F, A \f$
489// \return \f$ F^T A F\f$
490//
491template<typename T, Index N>
493Tensor<T, N>
495{
496 return t_dot(F, dot(A, F));
497}
498
499//
500// Push forward contravariant tensor
501// \param \f$ F, A \f$
502// \return \f$ F A F^T \f$
503//
504template<typename T, Index N>
506Tensor<T, N>
508{
509 return dot_t(dot(F, A), F);
510}
511
512//
513// Pull back contravariant tensor
514// \param \f$ F, A \f$
515// \return \f$ F^{-1} A F^{-T} \f$
516//
517template<typename T, Index N>
519Tensor<T, N>
521{
522 Index const
523 dimension = F.get_dimension();
524
526 G(dimension);
527
528 T const
529 J = det(F);
530
531 assert(J > 0.0);
532
533 switch (dimension) {
534
535 default:
536 MT_ERROR_EXIT("Supports only 2D and 3D.");
537 break;
538
539 case 3:
540 G(0,0) = (-F(1,2)*F(2,1) + F(1,1)*F(2,2)) / J;
541 G(0,1) = ( F(0,2)*F(2,1) - F(0,1)*F(2,2)) / J;
542 G(0,2) = (-F(0,2)*F(1,1) + F(0,1)*F(1,2)) / J;
543
544 G(1,0) = ( F(1,2)*F(2,0) - F(1,0)*F(2,2)) / J;
545 G(1,1) = (-F(0,2)*F(2,0) + F(0,0)*F(2,2)) / J;
546 G(1,2) = ( F(0,2)*F(1,0) - F(0,0)*F(1,2)) / J;
547
548 G(2,0) = (-F(1,1)*F(2,0) + F(1,0)*F(2,1)) / J;
549 G(2,1) = ( F(0,1)*F(2,0) - F(0,0)*F(2,1)) / J;
550 G(2,2) = (-F(0,1)*F(1,0) + F(0,0)*F(1,1)) / J;
551 break;
552
553 case 2:
554 G(0,0) = F(1,1) / J;
555 G(0,1) = -F(0,1) / J;
556
557 G(1,0) = -F(1,0) / J;
558 G(1,1) = F(0,0) / J;
559 break;
560
561 }
562
563 return dot_t(dot(G, A), G);
564}
565
566//
567// Piola transformation for vector
568// \param \f$ F, u \f$
569// \return \f$ \det F F^{-1} u \f$
570//
571template<typename T, Index N>
573Vector<T, N>
574piola(Tensor<T, N> const & F, Vector<T, N> const & u)
575{
576 Index const
577 dimension = F.get_dimension();
578
580 v(dimension);
581
582 switch (dimension) {
583
584 default:
585 MT_ERROR_EXIT("Supports only 2D and 3D.");
586 break;
587
588 case 3:
589 v(0) = (
590 (-F(1,2)*F(2,1) + F(1,1)*F(2,2)) * u(0) +
591 ( F(0,2)*F(2,1) - F(0,1)*F(2,2)) * u(1) +
592 (-F(0,2)*F(1,1) + F(0,1)*F(1,2)) * u(2));
593
594 v(1) = (
595 ( F(1,2)*F(2,0) - F(1,0)*F(2,2)) * u(0) +
596 (-F(0,2)*F(2,0) + F(0,0)*F(2,2)) * u(1) +
597 ( F(0,2)*F(1,0) - F(0,0)*F(1,2)) * u(2));
598
599 v(2) = (
600 (-F(1,1)*F(2,0) + F(1,0)*F(2,1)) * u(0) +
601 ( F(0,1)*F(2,0) - F(0,0)*F(2,1)) * u(1) +
602 (-F(0,1)*F(1,0) + F(0,0)*F(1,1)) * u(2));
603
604 break;
605
606 case 2:
607 v(0) = ( F(1,1) * u(0) - F(0,1) * u(1));
608 v(1) = (-F(1,0) * u(0) + F(0,0) * u(1));
609 break;
610
611 }
612
613 return v;
614}
615
616//
617// Inverse Piola transformation for vector
618// \param \f$ F, u \f$
619// \return \f$ (\det F)^{-1} F u \f$
620//
621template<typename T, Index N>
623Vector<T, N>
625{
626 Index const
627 dimension = F.get_dimension();
628
630 v(dimension);
631
632 T const
633 J = det(F);
634
635 assert(J > 0.0);
636
637 switch (dimension) {
638
639 default:
640 MT_ERROR_EXIT("Supports only 2D and 3D.");
641 break;
642
643 case 3:
644 v(0) = (F(0,0) * u(0) + F(0,1) * u(1) + F(0,2) * u(2)) / J;
645 v(1) = (F(1,0) * u(0) + F(1,1) * u(1) + F(1,2) * u(2)) / J;
646 v(2) = (F(2,0) * u(0) + F(2,1) * u(1) + F(2,2) * u(2)) / J;
647
648 break;
649
650 case 2:
651 v(0) = (F(0,0) * u(0) + F(0,1) * u(1)) / J;
652 v(1) = (F(1,0) * u(0) + F(1,1) * u(1)) / J;
653
654 break;
655
656 }
657
658 return v;
659}
660
661//
662// Piola transformation for tensor, applied on second
663// index. Useful for transforming Cauchy stress to 1PK stress.
664// \param \f$ F, \sigma \f$
665// \return \f$ \det F \sigma F^{-T} \f$
666//
667template<typename T, Index N>
669Tensor<T, N>
670piola(Tensor<T, N> const & F, Tensor<T, N> const & sigma)
671{
672 Index const
673 dimension = F.get_dimension();
674
676 G(dimension);
677
678 switch (dimension) {
679
680 default:
681 MT_ERROR_EXIT("Supports only 2D and 3D.");
682 break;
683
684 case 3:
685 G(0,0) = (-F(1,2)*F(2,1) + F(1,1)*F(2,2));
686 G(0,1) = ( F(0,2)*F(2,1) - F(0,1)*F(2,2));
687 G(0,2) = (-F(0,2)*F(1,1) + F(0,1)*F(1,2));
688
689 G(1,0) = ( F(1,2)*F(2,0) - F(1,0)*F(2,2));
690 G(1,1) = (-F(0,2)*F(2,0) + F(0,0)*F(2,2));
691 G(1,2) = ( F(0,2)*F(1,0) - F(0,0)*F(1,2));
692
693 G(2,0) = (-F(1,1)*F(2,0) + F(1,0)*F(2,1));
694 G(2,1) = ( F(0,1)*F(2,0) - F(0,0)*F(2,1));
695 G(2,2) = (-F(0,1)*F(1,0) + F(0,0)*F(1,1));
696 break;
697
698 case 2:
699 G(0,0) = F(1,1);
700 G(0,1) = -F(0,1);
701
702 G(1,0) = -F(1,0);
703 G(1,1) = F(0,0);
704 break;
705
706 }
707
708 return dot_t(sigma, G);
709}
710
711//
712// Inverse Piola transformation for tensor, applied on second
713// index. Useful for transforming 1PK stress to Cauchy stress.
714// \param \f$ F, P \f$
715// \return \f$ (\det F)^{-1} P F^T \f$
716//
717template<typename T, Index N>
719Tensor<T, N>
721{
722 T const
723 J = det(F);
724
725 assert(J > 0.0);
726
727 return dot_t(P, F) / J;
728}
729
730//
731// Smallest eigenvalue by inverse iteration.
732//
733template<typename T, Index N>
735T
737{
739 B = inverse(A);
740
741 T const
742 tolerance = machine_epsilon<T>();
743
744 Index const
745 dimension = A.get_dimension();
746
748 v(dimension, Filler::ONES);
749
750 Index const
751 maximum_iterations = 128;
752
753 T
754 relative_error = 1.0;
755
756 Index
757 k = 0;
758
759 while (relative_error > tolerance && k < maximum_iterations) {
760
761 Vector<T, N> const
762 w = v;
763
764 v = unit(B * w);
765
766 relative_error = norm(v - w) / norm(w);
767
768 ++k;
769 }
770
771 return v * A * v;
772}
773
774//
775// Check strict ellipticity condition for 4th-order tensor.
776// Assume A has major symmetries.
777//
778template<typename T, Index N>
780bool
782{
783 // Convert to 2nd-order tensor
785 B(A);
786
787 // Check bounds for eigenvalues
788 T const
789 lower_bound = bounds_eigenvalues(B).first;
790
791 if (lower_bound > 0.0) {
792 return true;
793 }
794
795 // Get eigenvalue closest to zero only
796 T const
797 smallest_eigenvalue = smallest_eigenvavlue(B);
798
799 if (smallest_eigenvalue > 0.0) {
800 return true;
801 }
802
803 return false;
804}
805
806//
807// Check strong ellipticity condition for 4th-order tensor.
808// Assume A has major and minor symmetries.
809//
810template<typename T, Index N>
811std::pair<bool, Vector<T, N>>
813{
814 bool
815 is_elliptic = true;
816
817 Index const
818 dimension = A.get_dimension();
819
821 eigenvector(dimension, Filler::ONES);
822
823 eigenvector /= dimension;
824
825 Index const
826 maximum_iterarions = 128;
827
828 T const
829 tolerance = machine_epsilon<T>();
830
831 T
832 error = 1.0;
833
834#if defined(KOKKOS_ENABLE_CUDA)
835 T
836 prev_eigenvalue = DBL_MAX;
837#else
838 using S = typename Sacado::ScalarType<T>::type;
839
840 T
841 prev_eigenvalue = std::numeric_limits<S>::max();
842#endif
843
844 T
845 curr_eigenvalue = prev_eigenvalue;
846
847 Index
848 iteration = 0;
849
850 while (error > tolerance && iteration < maximum_iterarions) {
851
853 Q = dot2(eigenvector, dot(A, eigenvector));
854
856 V;
857
859 D;
860
861 std::tie(V, D) = eig_sym(Q);
862
863 curr_eigenvalue = D(dimension - 1, dimension - 1);
864
865 eigenvector = col(V, dimension - 1);
866
867 error = std::abs(prev_eigenvalue) / std::abs(curr_eigenvalue) - 1.0;
868
869 prev_eigenvalue = curr_eigenvalue;
870
871 ++iteration;
872 }
873
874 if (curr_eigenvalue <= 0.0) {
875 is_elliptic = false;
876 }
877
878 return std::make_pair(is_elliptic, eigenvector);
879}
880
882} // namespace minitensor
883
884#endif // MiniTensor_Mechanics_h
#define KOKKOS_INLINE_FUNCTION
#define MT_ERROR_EXIT(...)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Vector< T, M > col(Matrix< T, M, N > const &A, Index const j)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > dot_t(Matrix< S, M, P > const &A, Matrix< T, N, P > const &B)
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > dot(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION Matrix< typename Promote< S, T >::type, M, N > t_dot(Matrix< S, P, M > const &A, Matrix< T, P, N > const &B)
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > dot2(Tensor3< T, N > const &A, Vector< S > const &u)
std::pair< Tensor< T, N >, Tensor< T, N > > eig_sym(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Tensor< T, N > inverse(Tensor< T, N > const &A)
std::pair< bool, Vector< T, N > > check_strong_ellipticity(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, N > piola(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION Vector< T, N > push_forward_contravariant(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION Tensor< T, N > dev(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, N > push_forward_covariant(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION Vector< T, N > pull_back_covariant(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION T smallest_eigenvalue(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Tensor< T, N > vol(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, N > pull_back_contravariant(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION Vector< T, N > piola_inverse(Tensor< T, N > const &F, Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION bool check_strict_ellipticity(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION T norm(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION T det(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION T trace(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Quaternion< T > unit(Quaternion< T > const &q)
uint32_t Index
Indexing type.