10#if !defined(MiniTensor_Mechanics_h)
11#define MiniTensor_Mechanics_h
25template<
typename T, Index N>
28vol(Tensor<T, N>
const & A);
35template<
typename T, Index N>
38dev(Tensor<T, N>
const & A);
46template<
typename T, Index N>
57template<
typename T, Index N>
68template<
typename T, Index N>
79template<
typename T, Index N>
90template<
typename T, Index N>
101template<
typename T, Index N>
112template<
typename T, Index N>
123template<
typename T, Index N>
134template<
typename T, Index N>
137piola(Tensor<T, N>
const & F, Vector<T, N>
const & u);
145template<
typename T, Index N>
148piola_inverse(Tensor<T, N>
const & F, Vector<T, N>
const & u);
157template<
typename T, Index N>
160piola(Tensor<T, N>
const & F, Tensor<T, N>
const & sigma);
169template<
typename T, Index N>
172piola_inverse(Tensor<T, N>
const & F, Tensor<T, N>
const & P);
177template<
typename T, Index N>
190template<
typename T, Index N>
201template<
typename T, Index N>
202std::pair<bool, Vector<T, N>>
213template<
typename T, Index N>
222 theta = (1.0/dimension) *
trace(A);
224 return theta * eye<T, N>(dimension);
231template<
typename T, Index N>
247template<
typename T, Index N>
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;
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;
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;
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;
302template<
typename T, Index N>
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);
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);
342template<
typename T, Index N>
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);
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);
382template<
typename T, Index N>
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;
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;
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;
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;
437template<
typename T, Index N>
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;
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;
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;
475 G(0,1) = -F(0,1) / J;
477 G(1,0) = -F(1,0) / J;
491template<
typename T, Index N>
504template<
typename T, Index N>
517template<
typename T, Index N>
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;
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;
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;
555 G(0,1) = -F(0,1) / J;
557 G(1,0) = -F(1,0) / J;
571template<
typename T, Index N>
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));
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));
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));
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));
621template<
typename T, Index N>
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;
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;
667template<
typename T, Index N>
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));
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));
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));
708 return dot_t(sigma, G);
717template<
typename T, Index N>
727 return dot_t(P, F) / J;
733template<
typename T, Index N>
742 tolerance = machine_epsilon<T>();
751 maximum_iterations = 128;
754 relative_error = 1.0;
759 while (relative_error > tolerance && k < maximum_iterations) {
766 relative_error =
norm(v - w) /
norm(w);
778template<
typename T, Index N>
789 lower_bound = bounds_eigenvalues(B).first;
791 if (lower_bound > 0.0) {
810template<
typename T, Index N>
811std::pair<bool, Vector<T, N>>
823 eigenvector /= dimension;
826 maximum_iterarions = 128;
829 tolerance = machine_epsilon<T>();
834#if defined(KOKKOS_ENABLE_CUDA)
836 prev_eigenvalue = DBL_MAX;
838 using S =
typename Sacado::ScalarType<T>::type;
841 prev_eigenvalue = std::numeric_limits<S>::max();
845 curr_eigenvalue = prev_eigenvalue;
850 while (error > tolerance && iteration < maximum_iterarions) {
853 Q =
dot2(eigenvector,
dot(A, eigenvector));
863 curr_eigenvalue = D(dimension - 1, dimension - 1);
865 eigenvector =
col(V, dimension - 1);
867 error = std::abs(prev_eigenvalue) / std::abs(curr_eigenvalue) - 1.0;
869 prev_eigenvalue = curr_eigenvalue;
874 if (curr_eigenvalue <= 0.0) {
878 return std::make_pair(is_elliptic, eigenvector);
#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.